Comparative analysis of Alignment/Mapping Pipelines in SARS-CoV-2-infected neuron cells
Mira Sohn
Raw data
- Study summary: Neurological complications are common in patients with COVID-19. While SARS-CoV-2, the causal pathogen of COVID-19, has been detected in some patient brains, its ability to infect brain cells and impact their function are not well understood, and experimental models using human brain cells are urgently needed. Here we investigated the susceptibility of human induced pluripotent stem cell (hiPSC)-derived monolayer brain cells and region-specific brain organoids to SARS-CoV-2 infection. We found modest numbers of infected neurons and astrocytes, but greater infection of choroid plexus epithelial cells. We optimized a protocol to generate choroid plexus organoids from hiPSCs, which revealed productive SARS-CoV-2 infection that leads to increased cell death and transcriptional dysregulation indicative of an inflammatory response and cellular function deficits. Together, our results provide evidence for SARS-CoV-2 neurotropism and support use of hiPSC-derived brain organoids as a platform to investigate the cellular susceptibility, disease mechanisms, and treatment strategies for SARS-CoV-2 infection. Bulk RNA-seq of choroid plexus organoids (CPOs) was performed on mock 72 hours post-infection (hpi), SARS-CoV-2 24 hpi, and SARS-CoV-2 72 hpi samples. All conditions were profiled in triplicate.
Loading packages
- AnnotationHub: https://bioconductor.org/packages/devel/bioc/vignettes/AnnotationHub/inst/doc/AnnotationHub.html
- DESeq2: http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
library(data.table)
library(tidyverse)
library(rmarkdown)
library(ggplot2)
library(pheatmap)
library(AnnotationHub)
library(tximport)
library(Rsubread)
library(DESeq2)
library(UpSetR)
library(ensembldb)Setting AnnotationHub
Assign your species of interest
DB <- "EnsDb" # Set your DB of interest
AnnotationSpecies <- "Homo sapiens" # Set your species
ah <- AnnotationHub(hub=getAnnotationHubOption("URL")) # Bring annotation DBRunning AnnotationHub
# Filter annotation of interest
ahQuery <- query(ah,
pattern=c(DB, AnnotationSpecies),
ignore.case=T)
# Select the most recent data
DBName <- mcols(ahQuery) %>%
rownames() %>%
tail(1)
AnnoDb <- ah[[DBName]]
# Explore your EnsDb object with following accessors:
# columns(AnnpDb)
# keytypes(AnnoDb)
# keys(AnnoDb, keytype=..)
# select(AnnoDb, keys=.., columns=.., keytype=...)
AnnoKey <- keys(AnnoDb, keytype="TXID")
# Note: Annotation has to be done with not genome but transcripts
AnnoDb <- select(AnnoDb,
AnnoKey,
keytype="TXID",
columns=c("GENEID", "GENENAME"))
# Check if your AnnoDb has been extracted and saved correctely
class(AnnoDb)## [1] "data.frame"head(AnnoDb) # The column 1 has to assign transcript (e.g. ENSEMBLTRANS)## TXID GENEID GENENAME
## 1 ENST00000000233 ENSG00000004059 ARF5
## 2 ENST00000000412 ENSG00000003056 M6PR
## 3 ENST00000000442 ENSG00000173153 ESRRA
## 4 ENST00000001008 ENSG00000004478 FKBP4
## 5 ENST00000001146 ENSG00000003137 CYP26B1
## 6 ENST00000002125 ENSG00000003509 NDUFAF7Metadata setting
# This code chunk needs to be written by yourself
# Define sample names
SampleNames <- c(paste0("Mock-rep", 1:3), paste0("SARS-CoV-2-rep", 1:3))
# Aligner names
Aligners <- c("Salmon", "STAR", "HISAT2")
# Define group level
GroupLevel <- c("Mock", "Covid")
# Define contrast for DE analysis
Contrast <- c("Group", GroupLevel)
# Set a function for file paths
path.fn <- function(head, tail) {
vec <- c(paste0("../",
head, # head = e.g. "hisat2", "star", or "salmon"
"_output/",
SampleNames,
tail)) # tail = file name after SampleNames
return(vec)
}
# Define .sf file path
sf <- path.fn("salmon",
".salmon_quant/quant.sf")
# Define STAR file path
star <- path.fn("star",
"Aligned.sortedByCoord.out.bam")
# Define HISAT2 file path
hisat <- path.fn("hisat2",
".sorted.bam")
# Define sample groups
group <- c(rep("Mock", 3), rep("Covid", 3))
# Create metadata
metadata <- data.frame(Sample=factor(SampleNames, levels=SampleNames),
Group=factor(group, levels=GroupLevel),
Salmon_path=sf,
STAR_path=star,
HISAT2_path=hisat)
# Assign row names with sample names
rownames(metadata) <- SampleNames
# Explore the metadata
print(metadata)## Sample Group
## Mock-rep1 Mock-rep1 Mock
## Mock-rep2 Mock-rep2 Mock
## Mock-rep3 Mock-rep3 Mock
## SARS-CoV-2-rep1 SARS-CoV-2-rep1 Covid
## SARS-CoV-2-rep2 SARS-CoV-2-rep2 Covid
## SARS-CoV-2-rep3 SARS-CoV-2-rep3 Covid
## Salmon_path
## Mock-rep1 ../salmon_output/Mock-rep1.salmon_quant/quant.sf
## Mock-rep2 ../salmon_output/Mock-rep2.salmon_quant/quant.sf
## Mock-rep3 ../salmon_output/Mock-rep3.salmon_quant/quant.sf
## SARS-CoV-2-rep1 ../salmon_output/SARS-CoV-2-rep1.salmon_quant/quant.sf
## SARS-CoV-2-rep2 ../salmon_output/SARS-CoV-2-rep2.salmon_quant/quant.sf
## SARS-CoV-2-rep3 ../salmon_output/SARS-CoV-2-rep3.salmon_quant/quant.sf
## STAR_path
## Mock-rep1 ../star_output/Mock-rep1Aligned.sortedByCoord.out.bam
## Mock-rep2 ../star_output/Mock-rep2Aligned.sortedByCoord.out.bam
## Mock-rep3 ../star_output/Mock-rep3Aligned.sortedByCoord.out.bam
## SARS-CoV-2-rep1 ../star_output/SARS-CoV-2-rep1Aligned.sortedByCoord.out.bam
## SARS-CoV-2-rep2 ../star_output/SARS-CoV-2-rep2Aligned.sortedByCoord.out.bam
## SARS-CoV-2-rep3 ../star_output/SARS-CoV-2-rep3Aligned.sortedByCoord.out.bam
## HISAT2_path
## Mock-rep1 ../hisat2_output/Mock-rep1.sorted.bam
## Mock-rep2 ../hisat2_output/Mock-rep2.sorted.bam
## Mock-rep3 ../hisat2_output/Mock-rep3.sorted.bam
## SARS-CoV-2-rep1 ../hisat2_output/SARS-CoV-2-rep1.sorted.bam
## SARS-CoV-2-rep2 ../hisat2_output/SARS-CoV-2-rep2.sorted.bam
## SARS-CoV-2-rep3 ../hisat2_output/SARS-CoV-2-rep3.sorted.bamfeatureCounts parameter setting
# "mm10", "mm9", "hg38", or "hg19"
annot.inbuilt <- "hg38"
# GTF file path
annot.ext <- "../../SEQC/reference_GENCODE/gencode.v36.primary_assembly.annotation.gtf"
# annotation type:
# e.g.: "gene_id", "transcript_id", or "gene_name"
GTF.attrType <- "gene_id"
# Number of cores
nthread <- 16
# Set a function importing counts from BAM files with featureCounts()
fcounts.fn <- function(vec) {
fc <- featureCounts(files=vec, # a vector assigning BAM file paths
annot.inbuilt=annot.inbuilt,
annot.ext=annot.ext,
GTF.attrType=GTF.attrType,
isGTFAnnotationFile=T,
nthread=nthread,
isPairedEnd=F, # Set this parameter correctly
verbose=T)
return(fc$counts)
}Importing counts
Importing Salmon counts
Note:
txi1 <- tximport(…, txOut=F)
txi2 <- tximport(…, txOut=T)
txi2 <- summarizedToGene(…)
counts extracted from txi1 and txi2 are the same
# Import gene level summarized counts
salmon.txi <- tximport(metadata$Salmon_path,
type = "salmon",
tx2gene=AnnoDb,
ignoreTxVersion=T,
txOut=F) # TRUE for transcript level, FALSE for gene level
# Extract the counts and save as a data frame
salmon.counts <- salmon.txi$counts
# Explore the salmon count data frame
head(salmon.counts)## [,1] [,2] [,3] [,4] [,5] [,6]
## ENSG00000000003 8991.942 7660.942 7504.000 11037.869 6193.989 7772.929
## ENSG00000000005 7.000 4.000 0.000 2.000 1.000 4.000
## ENSG00000000419 889.000 928.001 730.000 1211.001 777.000 991.000
## ENSG00000000457 699.297 564.374 566.745 832.783 437.987 596.733
## ENSG00000000460 452.703 366.157 397.262 740.218 388.013 514.268
## ENSG00000000938 0.000 1.000 0.000 0.000 0.000 0.000dim(salmon.counts)## [1] 60240 6summary(salmon.counts)## V1 V2 V3
## Min. : 0.0 Min. : 0.0 Min. : 0.0
## 1st Qu.: 0.0 1st Qu.: 0.0 1st Qu.: 0.0
## Median : 1.0 Median : 1.0 Median : 0.0
## Mean : 444.9 Mean : 434.5 Mean : 385.7
## 3rd Qu.: 21.0 3rd Qu.: 20.0 3rd Qu.: 16.0
## Max. :1809942.0 Max. :1898215.0 Max. :1628670.0
## V4 V5 V6
## Min. : 0.0 Min. : 0.0 Min. : 0.0
## 1st Qu.: 0.0 1st Qu.: 0.0 1st Qu.: 0.0
## Median : 1.0 Median : 0.0 Median : 0.0
## Mean : 562.4 Mean : 308.1 Mean : 427.5
## 3rd Qu.: 26.0 3rd Qu.: 14.0 3rd Qu.: 18.4
## Max. :1416833.0 Max. :775055.0 Max. :1265356.0Importing STAR counts
# Extract counts by running featureCounts()
star.counts <- fcounts.fn(metadata$STAR_path)##
## ========== _____ _ _ ____ _____ ______ _____
## ===== / ____| | | | _ \| __ \| ____| /\ | __ \
## ===== | (___ | | | | |_) | |__) | |__ / \ | | | |
## ==== \___ \| | | | _ <| _ /| __| / /\ \ | | | |
## ==== ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
## ========== |_____/ \____/|____/|_| \_\______/_/ \_\_____/
## Rsubread 2.4.0
##
## //========================== featureCounts setting ===========================\\
## || ||
## || Input files : 6 BAM files ||
## || ||
## || Mock-rep1Aligned.sortedByCoord.out.bam ||
## || Mock-rep2Aligned.sortedByCoord.out.bam ||
## || Mock-rep3Aligned.sortedByCoord.out.bam ||
## || SARS-CoV-2-rep1Aligned.sortedByCoord.out.bam ||
## || SARS-CoV-2-rep2Aligned.sortedByCoord.out.bam ||
## || SARS-CoV-2-rep3Aligned.sortedByCoord.out.bam ||
## || ||
## || Paired-end : no ||
## || Count read pairs : no ||
## || Annotation : gencode.v36.primary_assembly.annotation.gtf ... ||
## || Dir for temp files : . ||
## || Threads : 16 ||
## || Level : meta-feature level ||
## || Multimapping reads : counted ||
## || Multi-overlapping reads : not counted ||
## || Min overlapping bases : 1 ||
## || ||
## \\============================================================================//
##
## //================================= Running ==================================\\
## || ||
## || Load annotation file gencode.v36.primary_assembly.annotation.gtf ... ||
## || Features : 1430132 ||
## || Meta-features : 60719 ||
## || Chromosomes/contigs : 47 ||
## || ||
## || Process BAM file Mock-rep1Aligned.sortedByCoord.out.bam... ||
## || ||
## || Chromosomes/contigs in input file but not in annotation ||
## || KI270425.1 ||
## || KI270303.1 ||
## || KI270384.1 ||
## || KI270320.1 ||
## || KI270438.1 ||
## || KI270316.1 ||
## || KI270581.1 ||
## || KI270333.1 ||
## || KI270757.1 ||
## || KI270329.1 ||
## || KI270509.1 ||
## || KI270468.1 ||
## || KI270530.1 ||
## || KI270710.1 ||
## || KI270363.1 ||
## || KI270706.1 ||
## || KI270539.1 ||
## || KI270417.1 ||
## || KI270723.1 ||
## || KI270376.1 ||
## || KI270719.1 ||
## || KI270740.1 ||
## || KI270312.1 ||
## || KI270393.1 ||
## || KI270736.1 ||
## || KI270389.1 ||
## || GL000224.1 ||
## || KI270753.1 ||
## || KI270749.1 ||
## || KI270590.1 ||
## || KI270338.1 ||
## || KI270522.1 ||
## || KI270518.1 ||
## || KI270372.1 ||
## || KI270715.1 ||
## || KI270548.1 ||
## || KI270732.1 ||
## || KI270304.1 ||
## || KI270385.1 ||
## || KI270745.1 ||
## || KI270317.1 ||
## || KI270582.1 ||
## || KI270334.1 ||
## || KI270364.1 ||
## || KI270707.1 ||
## || KI270544.1 ||
## || KI270422.1 ||
## || KI270418.1 ||
## || KI270381.1 ||
## || KI270724.1 ||
## || KI270435.1 ||
## || GL000208.1 ||
## || KI270741.1 ||
## || KI270394.1 ||
## || KI270737.1 ||
## || KI270330.1 ||
## || KI270448.1 ||
## || KI270754.1 ||
## || KI270510.1 ||
## || KI270591.1 ||
## || KI270587.1 ||
## || KI270465.1 ||
## || KI270519.1 ||
## || KI270414.1 ||
## || KI270720.1 ||
## || KI270373.1 ||
## || KI270716.1 ||
## || KI270390.1 ||
## || KI270305.1 ||
## || KI270386.1 ||
## || KI270729.1 ||
## || GL000221.1 ||
## || KI270322.1 ||
## || KI270746.1 ||
## || KI270583.1 ||
## || KI270579.1 ||
## || KI270335.1 ||
## || KI270515.1 ||
## || KI270528.1 ||
## || KI270712.1 ||
## || KI270708.1 ||
## || KI270423.1 ||
## || KI270419.1 ||
## || KI270382.1 ||
## || KI270725.1 ||
## || KI270378.1 ||
## || KI270742.1 ||
## || KI270395.1 ||
## || KI270738.1 ||
## || GL000226.1 ||
## || KI270755.1 ||
## || KI270511.1 ||
## || KI270507.1 ||
## || KI270588.1 ||
## || KI270466.1 ||
## || GL000008.2 ||
## || KI270374.1 ||
## || KI270717.1 ||
## || KI270310.1 ||
## || KI270391.1 ||
## || KI270387.1 ||
## || KI270751.1 ||
## || KI270747.1 ||
## || KI270584.1 ||
## || KI270340.1 ||
## || KI270336.1 ||
## || KI270516.1 ||
## || KI270411.1 ||
## || KI270529.1 ||
## || KI270366.1 ||
## || KI270709.1 ||
## || KI270424.1 ||
## || KI270730.1 ||
## || KI270302.1 ||
## || KI270383.1 ||
## || KI270379.1 ||
## || GL000214.1 ||
## || KI270743.1 ||
## || KI270315.1 ||
## || KI270396.1 ||
## || KI270739.1 ||
## || KI270580.1 ||
## || KI270756.1 ||
## || KI270512.1 ||
## || KI270593.1 ||
## || KI270508.1 ||
## || KI270589.1 ||
## || KI270467.1 ||
## || KI270362.1 ||
## || KI270420.1 ||
## || KI270538.1 ||
## || KI270722.1 ||
## || KI270375.1 ||
## || KI270718.1 ||
## || KI270311.1 ||
## || KI270429.1 ||
## || KI270392.1 ||
## || KI270735.1 ||
## || KI270388.1 ||
## || KI270752.1 ||
## || KI270748.1 ||
## || KI270337.1 ||
## || KI270521.1 ||
## || KI270517.1 ||
## || KI270412.1 ||
## || KI270371.1 ||
## || KI270714.1 ||
## || ||
## || Single-end reads are included. ||
## || Total alignments : 54178083 ||
## || Successfully assigned alignments : 29102650 (53.7%) ||
## || Running time : 0.06 minutes ||
## || ||
## || Process BAM file Mock-rep2Aligned.sortedByCoord.out.bam... ||
## || ||
## || Chromosomes/contigs in input file but not in annotation ||
## || KI270425.1 ||
## || KI270303.1 ||
## || KI270384.1 ||
## || KI270320.1 ||
## || KI270438.1 ||
## || KI270316.1 ||
## || KI270581.1 ||
## || KI270333.1 ||
## || KI270757.1 ||
## || KI270329.1 ||
## || KI270509.1 ||
## || KI270468.1 ||
## || KI270530.1 ||
## || KI270710.1 ||
## || KI270363.1 ||
## || KI270706.1 ||
## || KI270539.1 ||
## || KI270417.1 ||
## || KI270723.1 ||
## || KI270376.1 ||
## || KI270719.1 ||
## || KI270740.1 ||
## || KI270312.1 ||
## || KI270393.1 ||
## || KI270736.1 ||
## || KI270389.1 ||
## || GL000224.1 ||
## || KI270753.1 ||
## || KI270749.1 ||
## || KI270590.1 ||
## || KI270338.1 ||
## || KI270522.1 ||
## || KI270518.1 ||
## || KI270372.1 ||
## || KI270715.1 ||
## || KI270548.1 ||
## || KI270732.1 ||
## || KI270304.1 ||
## || KI270385.1 ||
## || KI270745.1 ||
## || KI270317.1 ||
## || KI270582.1 ||
## || KI270334.1 ||
## || KI270364.1 ||
## || KI270707.1 ||
## || KI270544.1 ||
## || KI270422.1 ||
## || KI270418.1 ||
## || KI270381.1 ||
## || KI270724.1 ||
## || KI270435.1 ||
## || GL000208.1 ||
## || KI270741.1 ||
## || KI270394.1 ||
## || KI270737.1 ||
## || KI270330.1 ||
## || KI270448.1 ||
## || KI270754.1 ||
## || KI270510.1 ||
## || KI270591.1 ||
## || KI270587.1 ||
## || KI270465.1 ||
## || KI270519.1 ||
## || KI270414.1 ||
## || KI270720.1 ||
## || KI270373.1 ||
## || KI270716.1 ||
## || KI270390.1 ||
## || KI270305.1 ||
## || KI270386.1 ||
## || KI270729.1 ||
## || GL000221.1 ||
## || KI270322.1 ||
## || KI270746.1 ||
## || KI270583.1 ||
## || KI270579.1 ||
## || KI270335.1 ||
## || KI270515.1 ||
## || KI270528.1 ||
## || KI270712.1 ||
## || KI270708.1 ||
## || KI270423.1 ||
## || KI270419.1 ||
## || KI270382.1 ||
## || KI270725.1 ||
## || KI270378.1 ||
## || KI270742.1 ||
## || KI270395.1 ||
## || KI270738.1 ||
## || GL000226.1 ||
## || KI270755.1 ||
## || KI270511.1 ||
## || KI270507.1 ||
## || KI270588.1 ||
## || KI270466.1 ||
## || GL000008.2 ||
## || KI270374.1 ||
## || KI270717.1 ||
## || KI270310.1 ||
## || KI270391.1 ||
## || KI270387.1 ||
## || KI270751.1 ||
## || KI270747.1 ||
## || KI270584.1 ||
## || KI270340.1 ||
## || KI270336.1 ||
## || KI270516.1 ||
## || KI270411.1 ||
## || KI270529.1 ||
## || KI270366.1 ||
## || KI270709.1 ||
## || KI270424.1 ||
## || KI270730.1 ||
## || KI270302.1 ||
## || KI270383.1 ||
## || KI270379.1 ||
## || GL000214.1 ||
## || KI270743.1 ||
## || KI270315.1 ||
## || KI270396.1 ||
## || KI270739.1 ||
## || KI270580.1 ||
## || KI270756.1 ||
## || KI270512.1 ||
## || KI270593.1 ||
## || KI270508.1 ||
## || KI270589.1 ||
## || KI270467.1 ||
## || KI270362.1 ||
## || KI270420.1 ||
## || KI270538.1 ||
## || KI270722.1 ||
## || KI270375.1 ||
## || KI270718.1 ||
## || KI270311.1 ||
## || KI270429.1 ||
## || KI270392.1 ||
## || KI270735.1 ||
## || KI270388.1 ||
## || KI270752.1 ||
## || KI270748.1 ||
## || KI270337.1 ||
## || KI270521.1 ||
## || KI270517.1 ||
## || KI270412.1 ||
## || KI270371.1 ||
## || KI270714.1 ||
## || ||
## || Single-end reads are included. ||
## || Total alignments : 52009761 ||
## || Successfully assigned alignments : 28582595 (55.0%) ||
## || Running time : 0.06 minutes ||
## || ||
## || Process BAM file Mock-rep3Aligned.sortedByCoord.out.bam... ||
## || ||
## || Chromosomes/contigs in input file but not in annotation ||
## || KI270425.1 ||
## || KI270303.1 ||
## || KI270384.1 ||
## || KI270320.1 ||
## || KI270438.1 ||
## || KI270316.1 ||
## || KI270581.1 ||
## || KI270333.1 ||
## || KI270757.1 ||
## || KI270329.1 ||
## || KI270509.1 ||
## || KI270468.1 ||
## || KI270530.1 ||
## || KI270710.1 ||
## || KI270363.1 ||
## || KI270706.1 ||
## || KI270539.1 ||
## || KI270417.1 ||
## || KI270723.1 ||
## || KI270376.1 ||
## || KI270719.1 ||
## || KI270740.1 ||
## || KI270312.1 ||
## || KI270393.1 ||
## || KI270736.1 ||
## || KI270389.1 ||
## || GL000224.1 ||
## || KI270753.1 ||
## || KI270749.1 ||
## || KI270590.1 ||
## || KI270338.1 ||
## || KI270522.1 ||
## || KI270518.1 ||
## || KI270372.1 ||
## || KI270715.1 ||
## || KI270548.1 ||
## || KI270732.1 ||
## || KI270304.1 ||
## || KI270385.1 ||
## || KI270745.1 ||
## || KI270317.1 ||
## || KI270582.1 ||
## || KI270334.1 ||
## || KI270364.1 ||
## || KI270707.1 ||
## || KI270544.1 ||
## || KI270422.1 ||
## || KI270418.1 ||
## || KI270381.1 ||
## || KI270724.1 ||
## || KI270435.1 ||
## || GL000208.1 ||
## || KI270741.1 ||
## || KI270394.1 ||
## || KI270737.1 ||
## || KI270330.1 ||
## || KI270448.1 ||
## || KI270754.1 ||
## || KI270510.1 ||
## || KI270591.1 ||
## || KI270587.1 ||
## || KI270465.1 ||
## || KI270519.1 ||
## || KI270414.1 ||
## || KI270720.1 ||
## || KI270373.1 ||
## || KI270716.1 ||
## || KI270390.1 ||
## || KI270305.1 ||
## || KI270386.1 ||
## || KI270729.1 ||
## || GL000221.1 ||
## || KI270322.1 ||
## || KI270746.1 ||
## || KI270583.1 ||
## || KI270579.1 ||
## || KI270335.1 ||
## || KI270515.1 ||
## || KI270528.1 ||
## || KI270712.1 ||
## || KI270708.1 ||
## || KI270423.1 ||
## || KI270419.1 ||
## || KI270382.1 ||
## || KI270725.1 ||
## || KI270378.1 ||
## || KI270742.1 ||
## || KI270395.1 ||
## || KI270738.1 ||
## || GL000226.1 ||
## || KI270755.1 ||
## || KI270511.1 ||
## || KI270507.1 ||
## || KI270588.1 ||
## || KI270466.1 ||
## || GL000008.2 ||
## || KI270374.1 ||
## || KI270717.1 ||
## || KI270310.1 ||
## || KI270391.1 ||
## || KI270387.1 ||
## || KI270751.1 ||
## || KI270747.1 ||
## || KI270584.1 ||
## || KI270340.1 ||
## || KI270336.1 ||
## || KI270516.1 ||
## || KI270411.1 ||
## || KI270529.1 ||
## || KI270366.1 ||
## || KI270709.1 ||
## || KI270424.1 ||
## || KI270730.1 ||
## || KI270302.1 ||
## || KI270383.1 ||
## || KI270379.1 ||
## || GL000214.1 ||
## || KI270743.1 ||
## || KI270315.1 ||
## || KI270396.1 ||
## || KI270739.1 ||
## || KI270580.1 ||
## || KI270756.1 ||
## || KI270512.1 ||
## || KI270593.1 ||
## || KI270508.1 ||
## || KI270589.1 ||
## || KI270467.1 ||
## || KI270362.1 ||
## || KI270420.1 ||
## || KI270538.1 ||
## || KI270722.1 ||
## || KI270375.1 ||
## || KI270718.1 ||
## || KI270311.1 ||
## || KI270429.1 ||
## || KI270392.1 ||
## || KI270735.1 ||
## || KI270388.1 ||
## || KI270752.1 ||
## || KI270748.1 ||
## || KI270337.1 ||
## || KI270521.1 ||
## || KI270517.1 ||
## || KI270412.1 ||
## || KI270371.1 ||
## || KI270714.1 ||
## || ||
## || Single-end reads are included. ||
## || Total alignments : 44486311 ||
## || Successfully assigned alignments : 25165867 (56.6%) ||
## || Running time : 0.05 minutes ||
## || ||
## || Process BAM file SARS-CoV-2-rep1Aligned.sortedByCoord.out.bam... ||
## || ||
## || Chromosomes/contigs in input file but not in annotation ||
## || KI270425.1 ||
## || KI270303.1 ||
## || KI270384.1 ||
## || KI270320.1 ||
## || KI270438.1 ||
## || KI270316.1 ||
## || KI270581.1 ||
## || KI270333.1 ||
## || KI270757.1 ||
## || KI270329.1 ||
## || KI270509.1 ||
## || KI270468.1 ||
## || KI270530.1 ||
## || KI270710.1 ||
## || KI270363.1 ||
## || KI270706.1 ||
## || KI270539.1 ||
## || KI270417.1 ||
## || KI270723.1 ||
## || KI270376.1 ||
## || KI270719.1 ||
## || KI270740.1 ||
## || KI270312.1 ||
## || KI270393.1 ||
## || KI270736.1 ||
## || KI270389.1 ||
## || GL000224.1 ||
## || KI270753.1 ||
## || KI270749.1 ||
## || KI270590.1 ||
## || KI270338.1 ||
## || KI270522.1 ||
## || KI270518.1 ||
## || KI270372.1 ||
## || KI270715.1 ||
## || KI270548.1 ||
## || KI270732.1 ||
## || KI270304.1 ||
## || KI270385.1 ||
## || KI270745.1 ||
## || KI270317.1 ||
## || KI270582.1 ||
## || KI270334.1 ||
## || KI270364.1 ||
## || KI270707.1 ||
## || KI270544.1 ||
## || KI270422.1 ||
## || KI270418.1 ||
## || KI270381.1 ||
## || KI270724.1 ||
## || KI270435.1 ||
## || GL000208.1 ||
## || KI270741.1 ||
## || KI270394.1 ||
## || KI270737.1 ||
## || KI270330.1 ||
## || KI270448.1 ||
## || KI270754.1 ||
## || KI270510.1 ||
## || KI270591.1 ||
## || KI270587.1 ||
## || KI270465.1 ||
## || KI270519.1 ||
## || KI270414.1 ||
## || KI270720.1 ||
## || KI270373.1 ||
## || KI270716.1 ||
## || KI270390.1 ||
## || KI270305.1 ||
## || KI270386.1 ||
## || KI270729.1 ||
## || GL000221.1 ||
## || KI270322.1 ||
## || KI270746.1 ||
## || KI270583.1 ||
## || KI270579.1 ||
## || KI270335.1 ||
## || KI270515.1 ||
## || KI270528.1 ||
## || KI270712.1 ||
## || KI270708.1 ||
## || KI270423.1 ||
## || KI270419.1 ||
## || KI270382.1 ||
## || KI270725.1 ||
## || KI270378.1 ||
## || KI270742.1 ||
## || KI270395.1 ||
## || KI270738.1 ||
## || GL000226.1 ||
## || KI270755.1 ||
## || KI270511.1 ||
## || KI270507.1 ||
## || KI270588.1 ||
## || KI270466.1 ||
## || GL000008.2 ||
## || KI270374.1 ||
## || KI270717.1 ||
## || KI270310.1 ||
## || KI270391.1 ||
## || KI270387.1 ||
## || KI270751.1 ||
## || KI270747.1 ||
## || KI270584.1 ||
## || KI270340.1 ||
## || KI270336.1 ||
## || KI270516.1 ||
## || KI270411.1 ||
## || KI270529.1 ||
## || KI270366.1 ||
## || KI270709.1 ||
## || KI270424.1 ||
## || KI270730.1 ||
## || KI270302.1 ||
## || KI270383.1 ||
## || KI270379.1 ||
## || GL000214.1 ||
## || KI270743.1 ||
## || KI270315.1 ||
## || KI270396.1 ||
## || KI270739.1 ||
## || KI270580.1 ||
## || KI270756.1 ||
## || KI270512.1 ||
## || KI270593.1 ||
## || KI270508.1 ||
## || KI270589.1 ||
## || KI270467.1 ||
## || KI270362.1 ||
## || KI270420.1 ||
## || KI270538.1 ||
## || KI270722.1 ||
## || KI270375.1 ||
## || KI270718.1 ||
## || KI270311.1 ||
## || KI270429.1 ||
## || KI270392.1 ||
## || KI270735.1 ||
## || KI270388.1 ||
## || KI270752.1 ||
## || KI270748.1 ||
## || KI270337.1 ||
## || KI270521.1 ||
## || KI270517.1 ||
## || KI270412.1 ||
## || KI270371.1 ||
## || KI270714.1 ||
## || ||
## || Single-end reads are included. ||
## || Total alignments : 65798929 ||
## || Successfully assigned alignments : 35937122 (54.6%) ||
## || Running time : 0.07 minutes ||
## || ||
## || Process BAM file SARS-CoV-2-rep2Aligned.sortedByCoord.out.bam... ||
## || ||
## || Chromosomes/contigs in input file but not in annotation ||
## || KI270425.1 ||
## || KI270303.1 ||
## || KI270384.1 ||
## || KI270320.1 ||
## || KI270438.1 ||
## || KI270316.1 ||
## || KI270581.1 ||
## || KI270333.1 ||
## || KI270757.1 ||
## || KI270329.1 ||
## || KI270509.1 ||
## || KI270468.1 ||
## || KI270530.1 ||
## || KI270710.1 ||
## || KI270363.1 ||
## || KI270706.1 ||
## || KI270539.1 ||
## || KI270417.1 ||
## || KI270723.1 ||
## || KI270376.1 ||
## || KI270719.1 ||
## || KI270740.1 ||
## || KI270312.1 ||
## || KI270393.1 ||
## || KI270736.1 ||
## || KI270389.1 ||
## || GL000224.1 ||
## || KI270753.1 ||
## || KI270749.1 ||
## || KI270590.1 ||
## || KI270338.1 ||
## || KI270522.1 ||
## || KI270518.1 ||
## || KI270372.1 ||
## || KI270715.1 ||
## || KI270548.1 ||
## || KI270732.1 ||
## || KI270304.1 ||
## || KI270385.1 ||
## || KI270745.1 ||
## || KI270317.1 ||
## || KI270582.1 ||
## || KI270334.1 ||
## || KI270364.1 ||
## || KI270707.1 ||
## || KI270544.1 ||
## || KI270422.1 ||
## || KI270418.1 ||
## || KI270381.1 ||
## || KI270724.1 ||
## || KI270435.1 ||
## || GL000208.1 ||
## || KI270741.1 ||
## || KI270394.1 ||
## || KI270737.1 ||
## || KI270330.1 ||
## || KI270448.1 ||
## || KI270754.1 ||
## || KI270510.1 ||
## || KI270591.1 ||
## || KI270587.1 ||
## || KI270465.1 ||
## || KI270519.1 ||
## || KI270414.1 ||
## || KI270720.1 ||
## || KI270373.1 ||
## || KI270716.1 ||
## || KI270390.1 ||
## || KI270305.1 ||
## || KI270386.1 ||
## || KI270729.1 ||
## || GL000221.1 ||
## || KI270322.1 ||
## || KI270746.1 ||
## || KI270583.1 ||
## || KI270579.1 ||
## || KI270335.1 ||
## || KI270515.1 ||
## || KI270528.1 ||
## || KI270712.1 ||
## || KI270708.1 ||
## || KI270423.1 ||
## || KI270419.1 ||
## || KI270382.1 ||
## || KI270725.1 ||
## || KI270378.1 ||
## || KI270742.1 ||
## || KI270395.1 ||
## || KI270738.1 ||
## || GL000226.1 ||
## || KI270755.1 ||
## || KI270511.1 ||
## || KI270507.1 ||
## || KI270588.1 ||
## || KI270466.1 ||
## || GL000008.2 ||
## || KI270374.1 ||
## || KI270717.1 ||
## || KI270310.1 ||
## || KI270391.1 ||
## || KI270387.1 ||
## || KI270751.1 ||
## || KI270747.1 ||
## || KI270584.1 ||
## || KI270340.1 ||
## || KI270336.1 ||
## || KI270516.1 ||
## || KI270411.1 ||
## || KI270529.1 ||
## || KI270366.1 ||
## || KI270709.1 ||
## || KI270424.1 ||
## || KI270730.1 ||
## || KI270302.1 ||
## || KI270383.1 ||
## || KI270379.1 ||
## || GL000214.1 ||
## || KI270743.1 ||
## || KI270315.1 ||
## || KI270396.1 ||
## || KI270739.1 ||
## || KI270580.1 ||
## || KI270756.1 ||
## || KI270512.1 ||
## || KI270593.1 ||
## || KI270508.1 ||
## || KI270589.1 ||
## || KI270467.1 ||
## || KI270362.1 ||
## || KI270420.1 ||
## || KI270538.1 ||
## || KI270722.1 ||
## || KI270375.1 ||
## || KI270718.1 ||
## || KI270311.1 ||
## || KI270429.1 ||
## || KI270392.1 ||
## || KI270735.1 ||
## || KI270388.1 ||
## || KI270752.1 ||
## || KI270748.1 ||
## || KI270337.1 ||
## || KI270521.1 ||
## || KI270517.1 ||
## || KI270412.1 ||
## || KI270371.1 ||
## || KI270714.1 ||
## || ||
## || Single-end reads are included. ||
## || Total alignments : 34330831 ||
## || Successfully assigned alignments : 19826399 (57.8%) ||
## || Running time : 0.04 minutes ||
## || ||
## || Process BAM file SARS-CoV-2-rep3Aligned.sortedByCoord.out.bam... ||
## || ||
## || Chromosomes/contigs in input file but not in annotation ||
## || KI270425.1 ||
## || KI270303.1 ||
## || KI270384.1 ||
## || KI270320.1 ||
## || KI270438.1 ||
## || KI270316.1 ||
## || KI270581.1 ||
## || KI270333.1 ||
## || KI270757.1 ||
## || KI270329.1 ||
## || KI270509.1 ||
## || KI270468.1 ||
## || KI270530.1 ||
## || KI270710.1 ||
## || KI270363.1 ||
## || KI270706.1 ||
## || KI270539.1 ||
## || KI270417.1 ||
## || KI270723.1 ||
## || KI270376.1 ||
## || KI270719.1 ||
## || KI270740.1 ||
## || KI270312.1 ||
## || KI270393.1 ||
## || KI270736.1 ||
## || KI270389.1 ||
## || GL000224.1 ||
## || KI270753.1 ||
## || KI270749.1 ||
## || KI270590.1 ||
## || KI270338.1 ||
## || KI270522.1 ||
## || KI270518.1 ||
## || KI270372.1 ||
## || KI270715.1 ||
## || KI270548.1 ||
## || KI270732.1 ||
## || KI270304.1 ||
## || KI270385.1 ||
## || KI270745.1 ||
## || KI270317.1 ||
## || KI270582.1 ||
## || KI270334.1 ||
## || KI270364.1 ||
## || KI270707.1 ||
## || KI270544.1 ||
## || KI270422.1 ||
## || KI270418.1 ||
## || KI270381.1 ||
## || KI270724.1 ||
## || KI270435.1 ||
## || GL000208.1 ||
## || KI270741.1 ||
## || KI270394.1 ||
## || KI270737.1 ||
## || KI270330.1 ||
## || KI270448.1 ||
## || KI270754.1 ||
## || KI270510.1 ||
## || KI270591.1 ||
## || KI270587.1 ||
## || KI270465.1 ||
## || KI270519.1 ||
## || KI270414.1 ||
## || KI270720.1 ||
## || KI270373.1 ||
## || KI270716.1 ||
## || KI270390.1 ||
## || KI270305.1 ||
## || KI270386.1 ||
## || KI270729.1 ||
## || GL000221.1 ||
## || KI270322.1 ||
## || KI270746.1 ||
## || KI270583.1 ||
## || KI270579.1 ||
## || KI270335.1 ||
## || KI270515.1 ||
## || KI270528.1 ||
## || KI270712.1 ||
## || KI270708.1 ||
## || KI270423.1 ||
## || KI270419.1 ||
## || KI270382.1 ||
## || KI270725.1 ||
## || KI270378.1 ||
## || KI270742.1 ||
## || KI270395.1 ||
## || KI270738.1 ||
## || GL000226.1 ||
## || KI270755.1 ||
## || KI270511.1 ||
## || KI270507.1 ||
## || KI270588.1 ||
## || KI270466.1 ||
## || GL000008.2 ||
## || KI270374.1 ||
## || KI270717.1 ||
## || KI270310.1 ||
## || KI270391.1 ||
## || KI270387.1 ||
## || KI270751.1 ||
## || KI270747.1 ||
## || KI270584.1 ||
## || KI270340.1 ||
## || KI270336.1 ||
## || KI270516.1 ||
## || KI270411.1 ||
## || KI270529.1 ||
## || KI270366.1 ||
## || KI270709.1 ||
## || KI270424.1 ||
## || KI270730.1 ||
## || KI270302.1 ||
## || KI270383.1 ||
## || KI270379.1 ||
## || GL000214.1 ||
## || KI270743.1 ||
## || KI270315.1 ||
## || KI270396.1 ||
## || KI270739.1 ||
## || KI270580.1 ||
## || KI270756.1 ||
## || KI270512.1 ||
## || KI270593.1 ||
## || KI270508.1 ||
## || KI270589.1 ||
## || KI270467.1 ||
## || KI270362.1 ||
## || KI270420.1 ||
## || KI270538.1 ||
## || KI270722.1 ||
## || KI270375.1 ||
## || KI270718.1 ||
## || KI270311.1 ||
## || KI270429.1 ||
## || KI270392.1 ||
## || KI270735.1 ||
## || KI270388.1 ||
## || KI270752.1 ||
## || KI270748.1 ||
## || KI270337.1 ||
## || KI270521.1 ||
## || KI270517.1 ||
## || KI270412.1 ||
## || KI270371.1 ||
## || KI270714.1 ||
## || ||
## || Single-end reads are included. ||
## || Total alignments : 50803344 ||
## || Successfully assigned alignments : 27376989 (53.9%) ||
## || Running time : 0.06 minutes ||
## || ||
## || Write the final count table. ||
## || Write the read assignment summary. ||
## || ||
## \\============================================================================//# Explore the STAR count data frame
head(star.counts)## Mock-rep1Aligned.sortedByCoord.out.bam
## ENSG00000223972.5 24
## ENSG00000227232.5 52
## ENSG00000278267.1 0
## ENSG00000243485.5 8
## ENSG00000284332.1 0
## ENSG00000237613.2 327
## Mock-rep2Aligned.sortedByCoord.out.bam
## ENSG00000223972.5 13
## ENSG00000227232.5 46
## ENSG00000278267.1 0
## ENSG00000243485.5 9
## ENSG00000284332.1 0
## ENSG00000237613.2 311
## Mock-rep3Aligned.sortedByCoord.out.bam
## ENSG00000223972.5 17
## ENSG00000227232.5 32
## ENSG00000278267.1 0
## ENSG00000243485.5 7
## ENSG00000284332.1 0
## ENSG00000237613.2 270
## SARS-CoV-2-rep1Aligned.sortedByCoord.out.bam
## ENSG00000223972.5 29
## ENSG00000227232.5 39
## ENSG00000278267.1 0
## ENSG00000243485.5 5
## ENSG00000284332.1 0
## ENSG00000237613.2 101
## SARS-CoV-2-rep2Aligned.sortedByCoord.out.bam
## ENSG00000223972.5 6
## ENSG00000227232.5 26
## ENSG00000278267.1 0
## ENSG00000243485.5 2
## ENSG00000284332.1 0
## ENSG00000237613.2 97
## SARS-CoV-2-rep3Aligned.sortedByCoord.out.bam
## ENSG00000223972.5 16
## ENSG00000227232.5 45
## ENSG00000278267.1 0
## ENSG00000243485.5 10
## ENSG00000284332.1 0
## ENSG00000237613.2 116dim(star.counts)## [1] 60719 6summary(star.counts)## Mock-rep1Aligned.sortedByCoord.out.bam Mock-rep2Aligned.sortedByCoord.out.bam
## Min. : 0.0 Min. : 0.0
## 1st Qu.: 0.0 1st Qu.: 0.0
## Median : 2.0 Median : 1.0
## Mean : 479.3 Mean : 470.7
## 3rd Qu.: 31.0 3rd Qu.: 29.0
## Max. :2096551.0 Max. :2203523.0
## Mock-rep3Aligned.sortedByCoord.out.bam
## Min. : 0.0
## 1st Qu.: 0.0
## Median : 1.0
## Mean : 414.5
## 3rd Qu.: 23.0
## Max. :1850592.0
## SARS-CoV-2-rep1Aligned.sortedByCoord.out.bam
## Min. : 0.0
## 1st Qu.: 0.0
## Median : 2.0
## Mean : 591.9
## 3rd Qu.: 38.0
## Max. :1464055.0
## SARS-CoV-2-rep2Aligned.sortedByCoord.out.bam
## Min. : 0.0
## 1st Qu.: 0.0
## Median : 1.0
## Mean : 326.5
## 3rd Qu.: 20.0
## Max. :887660.0
## SARS-CoV-2-rep3Aligned.sortedByCoord.out.bam
## Min. : 0.0
## 1st Qu.: 0.0
## Median : 1.0
## Mean : 450.9
## 3rd Qu.: 27.0
## Max. :1410329.0Importing HISAT2 counts
# Extract counts by running featureCounts()
hisat2.counts <- fcounts.fn(metadata$HISAT2_path)##
## ========== _____ _ _ ____ _____ ______ _____
## ===== / ____| | | | _ \| __ \| ____| /\ | __ \
## ===== | (___ | | | | |_) | |__) | |__ / \ | | | |
## ==== \___ \| | | | _ <| _ /| __| / /\ \ | | | |
## ==== ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
## ========== |_____/ \____/|____/|_| \_\______/_/ \_\_____/
## Rsubread 2.4.0
##
## //========================== featureCounts setting ===========================\\
## || ||
## || Input files : 6 BAM files ||
## || ||
## || Mock-rep1.sorted.bam ||
## || Mock-rep2.sorted.bam ||
## || Mock-rep3.sorted.bam ||
## || SARS-CoV-2-rep1.sorted.bam ||
## || SARS-CoV-2-rep2.sorted.bam ||
## || SARS-CoV-2-rep3.sorted.bam ||
## || ||
## || Paired-end : no ||
## || Count read pairs : no ||
## || Annotation : gencode.v36.primary_assembly.annotation.gtf ... ||
## || Dir for temp files : . ||
## || Threads : 16 ||
## || Level : meta-feature level ||
## || Multimapping reads : counted ||
## || Multi-overlapping reads : not counted ||
## || Min overlapping bases : 1 ||
## || ||
## \\============================================================================//
##
## //================================= Running ==================================\\
## || ||
## || Load annotation file gencode.v36.primary_assembly.annotation.gtf ... ||
## || Features : 1430132 ||
## || Meta-features : 60719 ||
## || Chromosomes/contigs : 47 ||
## || ||
## || Process BAM file Mock-rep1.sorted.bam... ||
## || ||
## || Chromosomes/contigs in input file but not in annotation ||
## || KI270425.1 ||
## || KI270303.1 ||
## || KI270384.1 ||
## || KI270320.1 ||
## || KI270438.1 ||
## || KI270316.1 ||
## || KI270581.1 ||
## || KI270333.1 ||
## || KI270757.1 ||
## || KI270329.1 ||
## || KI270509.1 ||
## || KI270468.1 ||
## || KI270530.1 ||
## || KI270710.1 ||
## || KI270363.1 ||
## || KI270706.1 ||
## || KI270539.1 ||
## || KI270417.1 ||
## || KI270723.1 ||
## || KI270376.1 ||
## || KI270719.1 ||
## || KI270740.1 ||
## || KI270312.1 ||
## || KI270393.1 ||
## || KI270736.1 ||
## || KI270389.1 ||
## || GL000224.1 ||
## || KI270753.1 ||
## || KI270749.1 ||
## || KI270590.1 ||
## || KI270338.1 ||
## || KI270522.1 ||
## || KI270518.1 ||
## || KI270372.1 ||
## || KI270715.1 ||
## || KI270548.1 ||
## || KI270732.1 ||
## || KI270304.1 ||
## || KI270385.1 ||
## || KI270745.1 ||
## || KI270317.1 ||
## || KI270582.1 ||
## || KI270334.1 ||
## || KI270364.1 ||
## || KI270707.1 ||
## || KI270544.1 ||
## || KI270422.1 ||
## || KI270418.1 ||
## || KI270381.1 ||
## || KI270724.1 ||
## || KI270435.1 ||
## || GL000208.1 ||
## || KI270741.1 ||
## || KI270394.1 ||
## || KI270737.1 ||
## || KI270330.1 ||
## || KI270448.1 ||
## || KI270754.1 ||
## || KI270510.1 ||
## || KI270591.1 ||
## || KI270587.1 ||
## || KI270465.1 ||
## || KI270519.1 ||
## || KI270414.1 ||
## || KI270720.1 ||
## || KI270373.1 ||
## || KI270716.1 ||
## || KI270390.1 ||
## || KI270305.1 ||
## || KI270386.1 ||
## || KI270729.1 ||
## || GL000221.1 ||
## || KI270322.1 ||
## || KI270746.1 ||
## || KI270583.1 ||
## || KI270579.1 ||
## || KI270335.1 ||
## || KI270515.1 ||
## || KI270528.1 ||
## || KI270712.1 ||
## || KI270708.1 ||
## || KI270423.1 ||
## || KI270419.1 ||
## || KI270382.1 ||
## || KI270725.1 ||
## || KI270378.1 ||
## || KI270742.1 ||
## || KI270395.1 ||
## || KI270738.1 ||
## || GL000226.1 ||
## || KI270755.1 ||
## || KI270511.1 ||
## || KI270507.1 ||
## || KI270588.1 ||
## || KI270466.1 ||
## || GL000008.2 ||
## || KI270374.1 ||
## || KI270717.1 ||
## || KI270310.1 ||
## || KI270391.1 ||
## || KI270387.1 ||
## || KI270751.1 ||
## || KI270747.1 ||
## || KI270584.1 ||
## || KI270340.1 ||
## || KI270336.1 ||
## || KI270516.1 ||
## || KI270411.1 ||
## || KI270529.1 ||
## || KI270366.1 ||
## || KI270709.1 ||
## || KI270424.1 ||
## || KI270730.1 ||
## || KI270302.1 ||
## || KI270383.1 ||
## || KI270379.1 ||
## || GL000214.1 ||
## || KI270743.1 ||
## || KI270315.1 ||
## || KI270396.1 ||
## || KI270739.1 ||
## || KI270580.1 ||
## || KI270756.1 ||
## || KI270512.1 ||
## || KI270593.1 ||
## || KI270508.1 ||
## || KI270589.1 ||
## || KI270467.1 ||
## || KI270362.1 ||
## || KI270420.1 ||
## || KI270538.1 ||
## || KI270722.1 ||
## || KI270375.1 ||
## || KI270718.1 ||
## || KI270311.1 ||
## || KI270429.1 ||
## || KI270392.1 ||
## || KI270735.1 ||
## || KI270388.1 ||
## || KI270752.1 ||
## || KI270748.1 ||
## || KI270337.1 ||
## || KI270521.1 ||
## || KI270517.1 ||
## || KI270412.1 ||
## || KI270371.1 ||
## || KI270714.1 ||
## || ||
## || Single-end reads are included. ||
## || Total alignments : 54439819 ||
## || Successfully assigned alignments : 27242652 (50.0%) ||
## || Running time : 0.06 minutes ||
## || ||
## || Process BAM file Mock-rep2.sorted.bam... ||
## || ||
## || Chromosomes/contigs in input file but not in annotation ||
## || KI270425.1 ||
## || KI270303.1 ||
## || KI270384.1 ||
## || KI270320.1 ||
## || KI270438.1 ||
## || KI270316.1 ||
## || KI270581.1 ||
## || KI270333.1 ||
## || KI270757.1 ||
## || KI270329.1 ||
## || KI270509.1 ||
## || KI270468.1 ||
## || KI270530.1 ||
## || KI270710.1 ||
## || KI270363.1 ||
## || KI270706.1 ||
## || KI270539.1 ||
## || KI270417.1 ||
## || KI270723.1 ||
## || KI270376.1 ||
## || KI270719.1 ||
## || KI270740.1 ||
## || KI270312.1 ||
## || KI270393.1 ||
## || KI270736.1 ||
## || KI270389.1 ||
## || GL000224.1 ||
## || KI270753.1 ||
## || KI270749.1 ||
## || KI270590.1 ||
## || KI270338.1 ||
## || KI270522.1 ||
## || KI270518.1 ||
## || KI270372.1 ||
## || KI270715.1 ||
## || KI270548.1 ||
## || KI270732.1 ||
## || KI270304.1 ||
## || KI270385.1 ||
## || KI270745.1 ||
## || KI270317.1 ||
## || KI270582.1 ||
## || KI270334.1 ||
## || KI270364.1 ||
## || KI270707.1 ||
## || KI270544.1 ||
## || KI270422.1 ||
## || KI270418.1 ||
## || KI270381.1 ||
## || KI270724.1 ||
## || KI270435.1 ||
## || GL000208.1 ||
## || KI270741.1 ||
## || KI270394.1 ||
## || KI270737.1 ||
## || KI270330.1 ||
## || KI270448.1 ||
## || KI270754.1 ||
## || KI270510.1 ||
## || KI270591.1 ||
## || KI270587.1 ||
## || KI270465.1 ||
## || KI270519.1 ||
## || KI270414.1 ||
## || KI270720.1 ||
## || KI270373.1 ||
## || KI270716.1 ||
## || KI270390.1 ||
## || KI270305.1 ||
## || KI270386.1 ||
## || KI270729.1 ||
## || GL000221.1 ||
## || KI270322.1 ||
## || KI270746.1 ||
## || KI270583.1 ||
## || KI270579.1 ||
## || KI270335.1 ||
## || KI270515.1 ||
## || KI270528.1 ||
## || KI270712.1 ||
## || KI270708.1 ||
## || KI270423.1 ||
## || KI270419.1 ||
## || KI270382.1 ||
## || KI270725.1 ||
## || KI270378.1 ||
## || KI270742.1 ||
## || KI270395.1 ||
## || KI270738.1 ||
## || GL000226.1 ||
## || KI270755.1 ||
## || KI270511.1 ||
## || KI270507.1 ||
## || KI270588.1 ||
## || KI270466.1 ||
## || GL000008.2 ||
## || KI270374.1 ||
## || KI270717.1 ||
## || KI270310.1 ||
## || KI270391.1 ||
## || KI270387.1 ||
## || KI270751.1 ||
## || KI270747.1 ||
## || KI270584.1 ||
## || KI270340.1 ||
## || KI270336.1 ||
## || KI270516.1 ||
## || KI270411.1 ||
## || KI270529.1 ||
## || KI270366.1 ||
## || KI270709.1 ||
## || KI270424.1 ||
## || KI270730.1 ||
## || KI270302.1 ||
## || KI270383.1 ||
## || KI270379.1 ||
## || GL000214.1 ||
## || KI270743.1 ||
## || KI270315.1 ||
## || KI270396.1 ||
## || KI270739.1 ||
## || KI270580.1 ||
## || KI270756.1 ||
## || KI270512.1 ||
## || KI270593.1 ||
## || KI270508.1 ||
## || KI270589.1 ||
## || KI270467.1 ||
## || KI270362.1 ||
## || KI270420.1 ||
## || KI270538.1 ||
## || KI270722.1 ||
## || KI270375.1 ||
## || KI270718.1 ||
## || KI270311.1 ||
## || KI270429.1 ||
## || KI270392.1 ||
## || KI270735.1 ||
## || KI270388.1 ||
## || KI270752.1 ||
## || KI270748.1 ||
## || KI270337.1 ||
## || KI270521.1 ||
## || KI270517.1 ||
## || KI270412.1 ||
## || KI270371.1 ||
## || KI270714.1 ||
## || ||
## || Single-end reads are included. ||
## || Total alignments : 51357583 ||
## || Successfully assigned alignments : 26791035 (52.2%) ||
## || Running time : 0.06 minutes ||
## || ||
## || Process BAM file Mock-rep3.sorted.bam... ||
## || ||
## || Chromosomes/contigs in input file but not in annotation ||
## || KI270425.1 ||
## || KI270303.1 ||
## || KI270384.1 ||
## || KI270320.1 ||
## || KI270438.1 ||
## || KI270316.1 ||
## || KI270581.1 ||
## || KI270333.1 ||
## || KI270757.1 ||
## || KI270329.1 ||
## || KI270509.1 ||
## || KI270468.1 ||
## || KI270530.1 ||
## || KI270710.1 ||
## || KI270363.1 ||
## || KI270706.1 ||
## || KI270539.1 ||
## || KI270417.1 ||
## || KI270723.1 ||
## || KI270376.1 ||
## || KI270719.1 ||
## || KI270740.1 ||
## || KI270312.1 ||
## || KI270393.1 ||
## || KI270736.1 ||
## || KI270389.1 ||
## || GL000224.1 ||
## || KI270753.1 ||
## || KI270749.1 ||
## || KI270590.1 ||
## || KI270338.1 ||
## || KI270522.1 ||
## || KI270518.1 ||
## || KI270372.1 ||
## || KI270715.1 ||
## || KI270548.1 ||
## || KI270732.1 ||
## || KI270304.1 ||
## || KI270385.1 ||
## || KI270745.1 ||
## || KI270317.1 ||
## || KI270582.1 ||
## || KI270334.1 ||
## || KI270364.1 ||
## || KI270707.1 ||
## || KI270544.1 ||
## || KI270422.1 ||
## || KI270418.1 ||
## || KI270381.1 ||
## || KI270724.1 ||
## || KI270435.1 ||
## || GL000208.1 ||
## || KI270741.1 ||
## || KI270394.1 ||
## || KI270737.1 ||
## || KI270330.1 ||
## || KI270448.1 ||
## || KI270754.1 ||
## || KI270510.1 ||
## || KI270591.1 ||
## || KI270587.1 ||
## || KI270465.1 ||
## || KI270519.1 ||
## || KI270414.1 ||
## || KI270720.1 ||
## || KI270373.1 ||
## || KI270716.1 ||
## || KI270390.1 ||
## || KI270305.1 ||
## || KI270386.1 ||
## || KI270729.1 ||
## || GL000221.1 ||
## || KI270322.1 ||
## || KI270746.1 ||
## || KI270583.1 ||
## || KI270579.1 ||
## || KI270335.1 ||
## || KI270515.1 ||
## || KI270528.1 ||
## || KI270712.1 ||
## || KI270708.1 ||
## || KI270423.1 ||
## || KI270419.1 ||
## || KI270382.1 ||
## || KI270725.1 ||
## || KI270378.1 ||
## || KI270742.1 ||
## || KI270395.1 ||
## || KI270738.1 ||
## || GL000226.1 ||
## || KI270755.1 ||
## || KI270511.1 ||
## || KI270507.1 ||
## || KI270588.1 ||
## || KI270466.1 ||
## || GL000008.2 ||
## || KI270374.1 ||
## || KI270717.1 ||
## || KI270310.1 ||
## || KI270391.1 ||
## || KI270387.1 ||
## || KI270751.1 ||
## || KI270747.1 ||
## || KI270584.1 ||
## || KI270340.1 ||
## || KI270336.1 ||
## || KI270516.1 ||
## || KI270411.1 ||
## || KI270529.1 ||
## || KI270366.1 ||
## || KI270709.1 ||
## || KI270424.1 ||
## || KI270730.1 ||
## || KI270302.1 ||
## || KI270383.1 ||
## || KI270379.1 ||
## || GL000214.1 ||
## || KI270743.1 ||
## || KI270315.1 ||
## || KI270396.1 ||
## || KI270739.1 ||
## || KI270580.1 ||
## || KI270756.1 ||
## || KI270512.1 ||
## || KI270593.1 ||
## || KI270508.1 ||
## || KI270589.1 ||
## || KI270467.1 ||
## || KI270362.1 ||
## || KI270420.1 ||
## || KI270538.1 ||
## || KI270722.1 ||
## || KI270375.1 ||
## || KI270718.1 ||
## || KI270311.1 ||
## || KI270429.1 ||
## || KI270392.1 ||
## || KI270735.1 ||
## || KI270388.1 ||
## || KI270752.1 ||
## || KI270748.1 ||
## || KI270337.1 ||
## || KI270521.1 ||
## || KI270517.1 ||
## || KI270412.1 ||
## || KI270371.1 ||
## || KI270714.1 ||
## || ||
## || Single-end reads are included. ||
## || Total alignments : 44166286 ||
## || Successfully assigned alignments : 23859049 (54.0%) ||
## || Running time : 0.05 minutes ||
## || ||
## || Process BAM file SARS-CoV-2-rep1.sorted.bam... ||
## || ||
## || Chromosomes/contigs in input file but not in annotation ||
## || KI270425.1 ||
## || KI270303.1 ||
## || KI270384.1 ||
## || KI270320.1 ||
## || KI270438.1 ||
## || KI270316.1 ||
## || KI270581.1 ||
## || KI270333.1 ||
## || KI270757.1 ||
## || KI270329.1 ||
## || KI270509.1 ||
## || KI270468.1 ||
## || KI270530.1 ||
## || KI270710.1 ||
## || KI270363.1 ||
## || KI270706.1 ||
## || KI270539.1 ||
## || KI270417.1 ||
## || KI270723.1 ||
## || KI270376.1 ||
## || KI270719.1 ||
## || KI270740.1 ||
## || KI270312.1 ||
## || KI270393.1 ||
## || KI270736.1 ||
## || KI270389.1 ||
## || GL000224.1 ||
## || KI270753.1 ||
## || KI270749.1 ||
## || KI270590.1 ||
## || KI270338.1 ||
## || KI270522.1 ||
## || KI270518.1 ||
## || KI270372.1 ||
## || KI270715.1 ||
## || KI270548.1 ||
## || KI270732.1 ||
## || KI270304.1 ||
## || KI270385.1 ||
## || KI270745.1 ||
## || KI270317.1 ||
## || KI270582.1 ||
## || KI270334.1 ||
## || KI270364.1 ||
## || KI270707.1 ||
## || KI270544.1 ||
## || KI270422.1 ||
## || KI270418.1 ||
## || KI270381.1 ||
## || KI270724.1 ||
## || KI270435.1 ||
## || GL000208.1 ||
## || KI270741.1 ||
## || KI270394.1 ||
## || KI270737.1 ||
## || KI270330.1 ||
## || KI270448.1 ||
## || KI270754.1 ||
## || KI270510.1 ||
## || KI270591.1 ||
## || KI270587.1 ||
## || KI270465.1 ||
## || KI270519.1 ||
## || KI270414.1 ||
## || KI270720.1 ||
## || KI270373.1 ||
## || KI270716.1 ||
## || KI270390.1 ||
## || KI270305.1 ||
## || KI270386.1 ||
## || KI270729.1 ||
## || GL000221.1 ||
## || KI270322.1 ||
## || KI270746.1 ||
## || KI270583.1 ||
## || KI270579.1 ||
## || KI270335.1 ||
## || KI270515.1 ||
## || KI270528.1 ||
## || KI270712.1 ||
## || KI270708.1 ||
## || KI270423.1 ||
## || KI270419.1 ||
## || KI270382.1 ||
## || KI270725.1 ||
## || KI270378.1 ||
## || KI270742.1 ||
## || KI270395.1 ||
## || KI270738.1 ||
## || GL000226.1 ||
## || KI270755.1 ||
## || KI270511.1 ||
## || KI270507.1 ||
## || KI270588.1 ||
## || KI270466.1 ||
## || GL000008.2 ||
## || KI270374.1 ||
## || KI270717.1 ||
## || KI270310.1 ||
## || KI270391.1 ||
## || KI270387.1 ||
## || KI270751.1 ||
## || KI270747.1 ||
## || KI270584.1 ||
## || KI270340.1 ||
## || KI270336.1 ||
## || KI270516.1 ||
## || KI270411.1 ||
## || KI270529.1 ||
## || KI270366.1 ||
## || KI270709.1 ||
## || KI270424.1 ||
## || KI270730.1 ||
## || KI270302.1 ||
## || KI270383.1 ||
## || KI270379.1 ||
## || GL000214.1 ||
## || KI270743.1 ||
## || KI270315.1 ||
## || KI270396.1 ||
## || KI270739.1 ||
## || KI270580.1 ||
## || KI270756.1 ||
## || KI270512.1 ||
## || KI270593.1 ||
## || KI270508.1 ||
## || KI270589.1 ||
## || KI270467.1 ||
## || KI270362.1 ||
## || KI270420.1 ||
## || KI270538.1 ||
## || KI270722.1 ||
## || KI270375.1 ||
## || KI270718.1 ||
## || KI270311.1 ||
## || KI270429.1 ||
## || KI270392.1 ||
## || KI270735.1 ||
## || KI270388.1 ||
## || KI270752.1 ||
## || KI270748.1 ||
## || KI270337.1 ||
## || KI270521.1 ||
## || KI270517.1 ||
## || KI270412.1 ||
## || KI270371.1 ||
## || KI270714.1 ||
## || ||
## || Single-end reads are included. ||
## || Total alignments : 68677505 ||
## || Successfully assigned alignments : 33942278 (49.4%) ||
## || Running time : 0.08 minutes ||
## || ||
## || Process BAM file SARS-CoV-2-rep2.sorted.bam... ||
## || ||
## || Chromosomes/contigs in input file but not in annotation ||
## || KI270425.1 ||
## || KI270303.1 ||
## || KI270384.1 ||
## || KI270320.1 ||
## || KI270438.1 ||
## || KI270316.1 ||
## || KI270581.1 ||
## || KI270333.1 ||
## || KI270757.1 ||
## || KI270329.1 ||
## || KI270509.1 ||
## || KI270468.1 ||
## || KI270530.1 ||
## || KI270710.1 ||
## || KI270363.1 ||
## || KI270706.1 ||
## || KI270539.1 ||
## || KI270417.1 ||
## || KI270723.1 ||
## || KI270376.1 ||
## || KI270719.1 ||
## || KI270740.1 ||
## || KI270312.1 ||
## || KI270393.1 ||
## || KI270736.1 ||
## || KI270389.1 ||
## || GL000224.1 ||
## || KI270753.1 ||
## || KI270749.1 ||
## || KI270590.1 ||
## || KI270338.1 ||
## || KI270522.1 ||
## || KI270518.1 ||
## || KI270372.1 ||
## || KI270715.1 ||
## || KI270548.1 ||
## || KI270732.1 ||
## || KI270304.1 ||
## || KI270385.1 ||
## || KI270745.1 ||
## || KI270317.1 ||
## || KI270582.1 ||
## || KI270334.1 ||
## || KI270364.1 ||
## || KI270707.1 ||
## || KI270544.1 ||
## || KI270422.1 ||
## || KI270418.1 ||
## || KI270381.1 ||
## || KI270724.1 ||
## || KI270435.1 ||
## || GL000208.1 ||
## || KI270741.1 ||
## || KI270394.1 ||
## || KI270737.1 ||
## || KI270330.1 ||
## || KI270448.1 ||
## || KI270754.1 ||
## || KI270510.1 ||
## || KI270591.1 ||
## || KI270587.1 ||
## || KI270465.1 ||
## || KI270519.1 ||
## || KI270414.1 ||
## || KI270720.1 ||
## || KI270373.1 ||
## || KI270716.1 ||
## || KI270390.1 ||
## || KI270305.1 ||
## || KI270386.1 ||
## || KI270729.1 ||
## || GL000221.1 ||
## || KI270322.1 ||
## || KI270746.1 ||
## || KI270583.1 ||
## || KI270579.1 ||
## || KI270335.1 ||
## || KI270515.1 ||
## || KI270528.1 ||
## || KI270712.1 ||
## || KI270708.1 ||
## || KI270423.1 ||
## || KI270419.1 ||
## || KI270382.1 ||
## || KI270725.1 ||
## || KI270378.1 ||
## || KI270742.1 ||
## || KI270395.1 ||
## || KI270738.1 ||
## || GL000226.1 ||
## || KI270755.1 ||
## || KI270511.1 ||
## || KI270507.1 ||
## || KI270588.1 ||
## || KI270466.1 ||
## || GL000008.2 ||
## || KI270374.1 ||
## || KI270717.1 ||
## || KI270310.1 ||
## || KI270391.1 ||
## || KI270387.1 ||
## || KI270751.1 ||
## || KI270747.1 ||
## || KI270584.1 ||
## || KI270340.1 ||
## || KI270336.1 ||
## || KI270516.1 ||
## || KI270411.1 ||
## || KI270529.1 ||
## || KI270366.1 ||
## || KI270709.1 ||
## || KI270424.1 ||
## || KI270730.1 ||
## || KI270302.1 ||
## || KI270383.1 ||
## || KI270379.1 ||
## || GL000214.1 ||
## || KI270743.1 ||
## || KI270315.1 ||
## || KI270396.1 ||
## || KI270739.1 ||
## || KI270580.1 ||
## || KI270756.1 ||
## || KI270512.1 ||
## || KI270593.1 ||
## || KI270508.1 ||
## || KI270589.1 ||
## || KI270467.1 ||
## || KI270362.1 ||
## || KI270420.1 ||
## || KI270538.1 ||
## || KI270722.1 ||
## || KI270375.1 ||
## || KI270718.1 ||
## || KI270311.1 ||
## || KI270429.1 ||
## || KI270392.1 ||
## || KI270735.1 ||
## || KI270388.1 ||
## || KI270752.1 ||
## || KI270748.1 ||
## || KI270337.1 ||
## || KI270521.1 ||
## || KI270517.1 ||
## || KI270412.1 ||
## || KI270371.1 ||
## || KI270714.1 ||
## || ||
## || Single-end reads are included. ||
## || Total alignments : 35316473 ||
## || Successfully assigned alignments : 18758834 (53.1%) ||
## || Running time : 0.04 minutes ||
## || ||
## || Process BAM file SARS-CoV-2-rep3.sorted.bam... ||
## || ||
## || Chromosomes/contigs in input file but not in annotation ||
## || KI270425.1 ||
## || KI270303.1 ||
## || KI270384.1 ||
## || KI270320.1 ||
## || KI270438.1 ||
## || KI270316.1 ||
## || KI270581.1 ||
## || KI270333.1 ||
## || KI270757.1 ||
## || KI270329.1 ||
## || KI270509.1 ||
## || KI270468.1 ||
## || KI270530.1 ||
## || KI270710.1 ||
## || KI270363.1 ||
## || KI270706.1 ||
## || KI270539.1 ||
## || KI270417.1 ||
## || KI270723.1 ||
## || KI270376.1 ||
## || KI270719.1 ||
## || KI270740.1 ||
## || KI270312.1 ||
## || KI270393.1 ||
## || KI270736.1 ||
## || KI270389.1 ||
## || GL000224.1 ||
## || KI270753.1 ||
## || KI270749.1 ||
## || KI270590.1 ||
## || KI270338.1 ||
## || KI270522.1 ||
## || KI270518.1 ||
## || KI270372.1 ||
## || KI270715.1 ||
## || KI270548.1 ||
## || KI270732.1 ||
## || KI270304.1 ||
## || KI270385.1 ||
## || KI270745.1 ||
## || KI270317.1 ||
## || KI270582.1 ||
## || KI270334.1 ||
## || KI270364.1 ||
## || KI270707.1 ||
## || KI270544.1 ||
## || KI270422.1 ||
## || KI270418.1 ||
## || KI270381.1 ||
## || KI270724.1 ||
## || KI270435.1 ||
## || GL000208.1 ||
## || KI270741.1 ||
## || KI270394.1 ||
## || KI270737.1 ||
## || KI270330.1 ||
## || KI270448.1 ||
## || KI270754.1 ||
## || KI270510.1 ||
## || KI270591.1 ||
## || KI270587.1 ||
## || KI270465.1 ||
## || KI270519.1 ||
## || KI270414.1 ||
## || KI270720.1 ||
## || KI270373.1 ||
## || KI270716.1 ||
## || KI270390.1 ||
## || KI270305.1 ||
## || KI270386.1 ||
## || KI270729.1 ||
## || GL000221.1 ||
## || KI270322.1 ||
## || KI270746.1 ||
## || KI270583.1 ||
## || KI270579.1 ||
## || KI270335.1 ||
## || KI270515.1 ||
## || KI270528.1 ||
## || KI270712.1 ||
## || KI270708.1 ||
## || KI270423.1 ||
## || KI270419.1 ||
## || KI270382.1 ||
## || KI270725.1 ||
## || KI270378.1 ||
## || KI270742.1 ||
## || KI270395.1 ||
## || KI270738.1 ||
## || GL000226.1 ||
## || KI270755.1 ||
## || KI270511.1 ||
## || KI270507.1 ||
## || KI270588.1 ||
## || KI270466.1 ||
## || GL000008.2 ||
## || KI270374.1 ||
## || KI270717.1 ||
## || KI270310.1 ||
## || KI270391.1 ||
## || KI270387.1 ||
## || KI270751.1 ||
## || KI270747.1 ||
## || KI270584.1 ||
## || KI270340.1 ||
## || KI270336.1 ||
## || KI270516.1 ||
## || KI270411.1 ||
## || KI270529.1 ||
## || KI270366.1 ||
## || KI270709.1 ||
## || KI270424.1 ||
## || KI270730.1 ||
## || KI270302.1 ||
## || KI270383.1 ||
## || KI270379.1 ||
## || GL000214.1 ||
## || KI270743.1 ||
## || KI270315.1 ||
## || KI270396.1 ||
## || KI270739.1 ||
## || KI270580.1 ||
## || KI270756.1 ||
## || KI270512.1 ||
## || KI270593.1 ||
## || KI270508.1 ||
## || KI270589.1 ||
## || KI270467.1 ||
## || KI270362.1 ||
## || KI270420.1 ||
## || KI270538.1 ||
## || KI270722.1 ||
## || KI270375.1 ||
## || KI270718.1 ||
## || KI270311.1 ||
## || KI270429.1 ||
## || KI270392.1 ||
## || KI270735.1 ||
## || KI270388.1 ||
## || KI270752.1 ||
## || KI270748.1 ||
## || KI270337.1 ||
## || KI270521.1 ||
## || KI270517.1 ||
## || KI270412.1 ||
## || KI270371.1 ||
## || KI270714.1 ||
## || ||
## || Single-end reads are included. ||
## || Total alignments : 52224139 ||
## || Successfully assigned alignments : 26005814 (49.8%) ||
## || Running time : 0.06 minutes ||
## || ||
## || Write the final count table. ||
## || Write the read assignment summary. ||
## || ||
## \\============================================================================//# Explore the HISAT2 count data frame
head(hisat2.counts)## Mock-rep1.sorted.bam Mock-rep2.sorted.bam
## ENSG00000223972.5 20 11
## ENSG00000227232.5 46 62
## ENSG00000278267.1 0 0
## ENSG00000243485.5 6 8
## ENSG00000284332.1 0 0
## ENSG00000237613.2 267 256
## Mock-rep3.sorted.bam SARS-CoV-2-rep1.sorted.bam
## ENSG00000223972.5 12 21
## ENSG00000227232.5 36 57
## ENSG00000278267.1 0 0
## ENSG00000243485.5 4 2
## ENSG00000284332.1 0 0
## ENSG00000237613.2 231 87
## SARS-CoV-2-rep2.sorted.bam SARS-CoV-2-rep3.sorted.bam
## ENSG00000223972.5 4 10
## ENSG00000227232.5 35 57
## ENSG00000278267.1 0 0
## ENSG00000243485.5 2 6
## ENSG00000284332.1 0 0
## ENSG00000237613.2 81 89dim(hisat2.counts)## [1] 60719 6summary(hisat2.counts)## Mock-rep1.sorted.bam Mock-rep2.sorted.bam Mock-rep3.sorted.bam
## Min. : 0.0 Min. : 0.0 Min. : 0.0
## 1st Qu.: 0.0 1st Qu.: 0.0 1st Qu.: 0.0
## Median : 1.0 Median : 1.0 Median : 1.0
## Mean : 448.7 Mean : 441.2 Mean : 392.9
## 3rd Qu.: 27.0 3rd Qu.: 26.0 3rd Qu.: 21.0
## Max. :1876599.0 Max. :2004351.0 Max. :1678512.0
## SARS-CoV-2-rep1.sorted.bam SARS-CoV-2-rep2.sorted.bam
## Min. : 0 Min. : 0.0
## 1st Qu.: 0 1st Qu.: 0.0
## Median : 2 Median : 1.0
## Mean : 559 Mean : 308.9
## 3rd Qu.: 33 3rd Qu.: 17.0
## Max. :1410876 Max. :804577.0
## SARS-CoV-2-rep3.sorted.bam
## Min. : 0.0
## 1st Qu.: 0.0
## Median : 1.0
## Mean : 428.3
## 3rd Qu.: 24.0
## Max. :1328490.0Data cleaning: sample and gene annotation
countList <- list(salmon.counts,
star.counts,
hisat2.counts)
# Assign names of the count data frames in the count list
names(countList) <- Aligners
# Set a function cleaning the count data frame
clean.fn <- function(df) {
# Convert to a data frame
df <- as.data.frame(df)
# Assign column names
names(df) <- SampleNames
# Bring row names to a column
df <- df %>% rownames_to_column(var="GENEID")
return(df)
}
# Set a function to drop GENEID version
clean.annotation.fn <- function(df) {
# Re-annotate without version specification
df <- separate(df, "GENEID", c("GENEID", "Version"))
# Remove version column
df <- df[, colnames(df) != "Version"]
return(df)
}
# Move GENEID to a column
for (x in Aligners) {
countList[[x]] <- clean.fn(countList[[x]])
}
# Remove version of GENEID and duplicated rows in STAR & HISAT2 count tables
for (x in Aligners) {
countList[[x]] <- clean.annotation.fn(countList[[x]]) %>%
distinct()
}
# Explore the cleaned count data frames
head(countList[[1]])## GENEID Mock-rep1 Mock-rep2 Mock-rep3 SARS-CoV-2-rep1 SARS-CoV-2-rep2
## 1 ENSG00000000003 8991.942 7660.942 7504.000 11037.869 6193.989
## 2 ENSG00000000005 7.000 4.000 0.000 2.000 1.000
## 3 ENSG00000000419 889.000 928.001 730.000 1211.001 777.000
## 4 ENSG00000000457 699.297 564.374 566.745 832.783 437.987
## 5 ENSG00000000460 452.703 366.157 397.262 740.218 388.013
## 6 ENSG00000000938 0.000 1.000 0.000 0.000 0.000
## SARS-CoV-2-rep3
## 1 7772.929
## 2 4.000
## 3 991.000
## 4 596.733
## 5 514.268
## 6 0.000head(countList[[2]])## GENEID Mock-rep1 Mock-rep2 Mock-rep3 SARS-CoV-2-rep1 SARS-CoV-2-rep2
## 1 ENSG00000223972 24 13 17 29 6
## 2 ENSG00000227232 52 46 32 39 26
## 3 ENSG00000278267 0 0 0 0 0
## 4 ENSG00000243485 8 9 7 5 2
## 5 ENSG00000284332 0 0 0 0 0
## 6 ENSG00000237613 327 311 270 101 97
## SARS-CoV-2-rep3
## 1 16
## 2 45
## 3 0
## 4 10
## 5 0
## 6 116head(countList[[3]])## GENEID Mock-rep1 Mock-rep2 Mock-rep3 SARS-CoV-2-rep1 SARS-CoV-2-rep2
## 1 ENSG00000223972 20 11 12 21 4
## 2 ENSG00000227232 46 62 36 57 35
## 3 ENSG00000278267 0 0 0 0 0
## 4 ENSG00000243485 6 8 4 2 2
## 5 ENSG00000284332 0 0 0 0 0
## 6 ENSG00000237613 267 256 231 87 81
## SARS-CoV-2-rep3
## 1 10
## 2 57
## 3 0
## 4 6
## 5 0
## 6 89dim(countList[[1]])## [1] 60240 7dim(countList[[2]])## [1] 60675 7dim(countList[[3]])## [1] 60695 7sum(duplicated(countList[[1]]))## [1] 0sum(duplicated(countList[[2]]))## [1] 0sum(duplicated(countList[[3]]))## [1] 0# Convert Salmon counts to integers
countList[["Salmon"]] <- cbind(GENEID=countList[["Salmon"]][, "GENEID"],
round(countList[["Salmon"]][,
colnames(countList[["Salmon"]]) %in% SampleNames]))
# Explore the cleaned count data frames
head(countList[[1]])## GENEID Mock-rep1 Mock-rep2 Mock-rep3 SARS-CoV-2-rep1 SARS-CoV-2-rep2
## 1 ENSG00000000003 8992 7661 7504 11038 6194
## 2 ENSG00000000005 7 4 0 2 1
## 3 ENSG00000000419 889 928 730 1211 777
## 4 ENSG00000000457 699 564 567 833 438
## 5 ENSG00000000460 453 366 397 740 388
## 6 ENSG00000000938 0 1 0 0 0
## SARS-CoV-2-rep3
## 1 7773
## 2 4
## 3 991
## 4 597
## 5 514
## 6 0Plotting sequencing depth
Number of total counts per sample
# Set a function generating a data frame with sequencing depth
seq.depth.fn <- function(df, aligner) {
seqdf <- as.data.frame(colSums(df[, SampleNames])) %>%
rownames_to_column (var="Sample") %>%
mutate(Aligner=aligner)
names(seqdf) <- c("Sample", "Count", "Aligner")
return(seqdf)
}
# Set a function for a bar plot comparing values
comparing.barplot.fn <- function(df, yval, title, ytitle) {
ggplot(df,
aes(x=Sample, y=yval, group=Aligner, fill=Aligner)) +
geom_bar(stat="identity", position="dodge") +
theme_bw() +
ggtitle(title) +
ylab(ytitle)
}
# Initialize the seq depth data frame with the first aligner
seq.depth.df <- seq.depth.fn(countList[[1]], Aligners[1])
# Extend the seq depth data frame with the rest of aligners
for (x in Aligners) {
if (x %in% Aligners[2:length(Aligners)]) {
seq.depth.df <- rbind(seq.depth.df,
seq.depth.fn(countList[[x]], x))
}
}
# Explore how the data frame
print(seq.depth.df)## Sample Count Aligner
## 1 Mock-rep1 26801012 Salmon
## 2 Mock-rep2 26175554 Salmon
## 3 Mock-rep3 23236794 Salmon
## 4 SARS-CoV-2-rep1 33877808 Salmon
## 5 SARS-CoV-2-rep2 18560926 Salmon
## 6 SARS-CoV-2-rep3 25751215 Salmon
## 7 Mock-rep1 29098819 STAR
## 8 Mock-rep2 28578955 STAR
## 9 Mock-rep3 25163540 STAR
## 10 SARS-CoV-2-rep1 35933153 STAR
## 11 SARS-CoV-2-rep2 19824100 STAR
## 12 SARS-CoV-2-rep3 27373805 STAR
## 13 Mock-rep1 27242642 HISAT2
## 14 Mock-rep2 26791029 HISAT2
## 15 Mock-rep3 23859044 HISAT2
## 16 SARS-CoV-2-rep1 33942259 HISAT2
## 17 SARS-CoV-2-rep2 18758823 HISAT2
## 18 SARS-CoV-2-rep3 26005791 HISAT2summary(seq.depth.df)## Sample Count Aligner
## Length:18 Min. :18560926 Length:18
## Class :character 1st Qu.:24185168 Class :character
## Mode :character Median :26483292 Mode :character
## Mean :26498626
## 3rd Qu.:28277668
## Max. :35933153# Convert character vectors to factors
seq.depth.df$Sample <- factor(seq.depth.df$Sample,
levels=SampleNames)
seq.depth.df$Aligner <- factor(seq.depth.df$Aligner,
levels=Aligners)
# Create a plot presenting sequencing depth
comparing.barplot.fn(seq.depth.df,
seq.depth.df$Count,
"Sequencing Depth by Sample and Aligner",
"Count")
Generating DESeq2 objects
- vst() was run to perform variance stabilizing transformation instead of rlog() which takes longer time with similar characteristics.
- The vsd object created by vst() is used for not DE analysis but QC.
- References: DESeq2 doc “Transcript abundance files”, DESeq2 doc “Variance stabilizing transformation”
# Initialize new lists for storing dds objects
ddsList <- countList
# Initialize new lists for storing vsd objects
vsdList <- countList
for (x in Aligners) {
# Create a count matrix from the count data frame
m <- countList[[x]][, colnames(countList[[x]]) != "GENEID"] %>%
as.matrix()
# Assigne row names
rownames(m) <- countList[[x]]$GENEID
# Generate a DESeq2 object
ddsList[[x]] <- DESeqDataSetFromMatrix(m,
colData=metadata,
design=~Group)
# Conduct vst
vsdList[[x]] <- varianceStabilizingTransformation(ddsList[[x]],
blind=TRUE)
}
# Explore generated objects
summary(ddsList)## Length Class Mode
## Salmon 60240 DESeqDataSet S4
## STAR 60675 DESeqDataSet S4
## HISAT2 60695 DESeqDataSet S4summary(vsdList)## Length Class Mode
## Salmon 60240 DESeqTransform S4
## STAR 60675 DESeqTransform S4
## HISAT2 60695 DESeqTransform S4Estimating size factors
- black dashed line: size factor = 1
- Reference: DESeq2 doc “Sample-/gene-dependent normalization factors”
# Calculate and add size factors to the DEseq object
for (x in Aligners) {
ddsList[[x]] <- estimateSizeFactors(ddsList[[x]])
}
# Set a function summarizing size factors by aligner to a data frame
sfactor.fn <- function(df, aligner) {
sizefactor <- as.data.frame(round(sizeFactors(df), 3)) %>%
rownames_to_column(var="Sample") %>%
mutate(Aligner=aligner)
names(sizefactor) <- c("Sample", "Size_Factor", "Aligner")
return(sizefactor)
}
# Initialize a data frame with the first aligner
size.factor.df <- sfactor.fn(ddsList[[1]], Aligners[1])
for (x in Aligners) {
if (x != Aligners[1]) {
size.factor.df <- rbind(size.factor.df,
sfactor.fn(ddsList[[x]], x))
}
}
# Explore the data frame
print(size.factor.df)## Sample Size_Factor Aligner
## 1 Mock-rep1 1.095 Salmon
## 2 Mock-rep2 1.050 Salmon
## 3 Mock-rep3 0.889 Salmon
## 4 SARS-CoV-2-rep1 1.362 Salmon
## 5 SARS-CoV-2-rep2 0.754 Salmon
## 6 SARS-CoV-2-rep3 1.023 Salmon
## 7 Mock-rep1 1.101 STAR
## 8 Mock-rep2 1.056 STAR
## 9 Mock-rep3 0.886 STAR
## 10 SARS-CoV-2-rep1 1.363 STAR
## 11 SARS-CoV-2-rep2 0.747 STAR
## 12 SARS-CoV-2-rep3 1.020 STAR
## 13 Mock-rep1 1.091 HISAT2
## 14 Mock-rep2 1.048 HISAT2
## 15 Mock-rep3 0.894 HISAT2
## 16 SARS-CoV-2-rep1 1.359 HISAT2
## 17 SARS-CoV-2-rep2 0.748 HISAT2
## 18 SARS-CoV-2-rep3 1.025 HISAT2# Convert character vectors to factors
size.factor.df$Sample <- factor(size.factor.df$Sample,
levels=SampleNames)
size.factor.df$Aligner <- factor(size.factor.df$Aligner,
levels=Aligners)
# Plot calculated size factors
comparing.barplot.fn(size.factor.df,
size.factor.df$Size_Factor,
"Size Factors by Aligner and Sample",
"Size Factor") + geom_hline(yintercept=1, linetype="dashed", color="black", size=1)
Estimating dispersion and Wald test
- Dispersion is calculated as a measure of variation instead of variance since variance gets larger when gene expression gets higher.
- Wald test is the default setting of DESeq2 which tests null hypothesis between two groups. You should use Likelihood ratio test (LRT) when comparing more than two groups.
- References: Harvard Chan Bioinformatics Core workshop I, Harvard Chan Bioinformatics Core workshop II, Harvard Chan Bioinformatics Core workshop III, DESeq2 "Wald test indivisual steps, DESeq2 doc “Likelihood ratio test”
for (x in Aligners) {
# Dispersion
ddsList[[x]] <- estimateDispersions(ddsList[[x]])
# Wald test
ddsList[[x]] <- nbinomWaldTest(ddsList[[x]])
}
# Explore generated data in the dds object
ddsList[[1]]## class: DESeqDataSet
## dim: 60240 6
## metadata(1): version
## assays(4): counts mu H cooks
## rownames(60240): ENSG00000000003 ENSG00000000005 ... ENSG00000288674
## ENSG00000288675
## rowData names(22): baseMean baseVar ... deviance maxCooks
## colnames(6): Mock-rep1 Mock-rep2 ... SARS-CoV-2-rep2 SARS-CoV-2-rep3
## colData names(6): Sample Group ... HISAT2_path sizeFactorSample QC: Principal Component Analysis (PCA)
Identifies source of variation and sample outliers
- Reference: DESeq2 doc “Principal component plot of the samples”, DESeq2 doc “Heatmap of the sample-to-sample distances”
# Assigne what to compare
GroupOfInterest <- Contrast[1]
# Set a function for sample pca
qcpca.fn <- function(obj, title) {
plotPCA(obj,
intgroup=GroupOfInterest,
returnData=FALSE) + theme_bw() + ggtitle(paste("PCA:", title))
}
# Print the plots
qcpca.fn(vsdList[[1]], Aligners[1]) 
qcpca.fn(vsdList[[2]], Aligners[2])
qcpca.fn(vsdList[[3]], Aligners[3]) 
Sample QC: Sample Correlation Heatmap
Identifies distance between samples & correlation in a group
# Heatmap annotation
HeatmapAnno <- metadata[, c("Sample", "Group")]
# Set a function generating a correlation heatmap
cheatmap.fn <- function(df, title) {
# Extract a normalized count matrix
vm <- assay(df)
# Generate a correlation matrix
cm <- cor(vm)
# Generate a heatmap
pheatmap(cm,
annotation=HeatmapAnno,
main=paste("Sample Correlation Heatmap:", title))
}
# Print the heatmaps
cheatmap.fn(vsdList[[1]], Aligners[1])
cheatmap.fn(vsdList[[2]], Aligners[2])
cheatmap.fn(vsdList[[3]], Aligners[3])
Running DE analysis
# Run DESeq
for (x in Aligners) {
ddsList[[x]] <- DESeq(ddsList[[x]])
# Check result names
ResNames <- resultsNames(ddsList[[x]])
print(ResNames)
}## [1] "Intercept" "Group_Covid_vs_Mock"
## [1] "Intercept" "Group_Covid_vs_Mock"
## [1] "Intercept" "Group_Covid_vs_Mock"Creating dispersion plots
- Dispersion is important since estimation by DESeq2 algorithm is based on the assumption that genes with similar expression levels have similar dispersion. If an RNA-seq dataset doesn’t satisfy this assumption, use other DE algorithms than DESeq2.
- References: DESeq2 doc "Dispersion plot and fitting alternatives, Harvard Chan Bioinformatics Core workshop
# Set a function plotting dispersion
dplot.fn <- function(dds, title) {
plotDispEsts(dds,
main=paste("Dispersion over Counts:", title))
}
# Plot dispersion patterns
dplot.fn(ddsList[[1]], Aligners[1])
dplot.fn(ddsList[[2]], Aligners[2])
dplot.fn(ddsList[[3]], Aligners[3])
# Do they fit well with the DESeq2 estimation model?Setting how to extract fold-change results
Change variables below
- The alpha denotes threshold of false discovery rate (FDR) assigned by users.
- In this analysis, the alpha is set to 0.1
# Set FDR threshold alpha
alpha=0.1
# Set the coefficients to compare
Coef <- ResNames[-1]
print(Coef) ## [1] "Group_Covid_vs_Mock"# Set a function to clean a result table
lfctable.fn <- function(df) {
df <- df %>%
rownames_to_column(var="GENEID") %>%
mutate(FDR=ifelse(padj < 0.1 & !is.na(padj),
"< 0.1",
"> 0.1"))
return(df)
}
# Set a function extracting results
extract.lfc.fn <- function(dds) {
res <- results(dds, contrast=Contrast, alpha=alpha)
lfctable.fn(as.data.frame(res))
return(lfctable.fn(as.data.frame(res)))
}Extracting log2FoldChanges
You can change alpha depending on your interest of FDR level
Shrinkage is NOT applied in this analysis
# Initialize a list storing lfc data frames
lfcList <- countList
# Extract DE results
# The Contrast variable was defined in the previous chunk
# Extraction with no shrinkage
# alpha: FDR threshold
for (x in Aligners) {
lfcList[[x]] <- extract.lfc.fn(ddsList[[x]]) %>% mutate(Alignment=x)
print(head(lfcList[[x]]))
}## GENEID baseMean log2FoldChange lfcSE stat pvalue
## 1 ENSG00000000003 7978.760054 0.001711158 0.1039000 0.01646927 0.98686002
## 2 ENSG00000000005 2.817915 0.617944320 1.5219706 0.40601594 0.68473089
## 3 ENSG00000000419 900.991910 -0.196869661 0.1265716 -1.55540134 0.11985051
## 4 ENSG00000000457 598.304194 0.028506669 0.1403096 0.20316984 0.83900228
## 5 ENSG00000000460 461.594729 -0.370406019 0.1542869 -2.40076109 0.01636101
## 6 ENSG00000000938 0.158719 0.974794073 4.0804729 0.23889243 0.81118900
## padj FDR Alignment
## 1 0.99521804 > 0.1 Salmon
## 2 NA > 0.1 Salmon
## 3 0.33075557 > 0.1 Salmon
## 4 0.93997529 > 0.1 Salmon
## 5 0.08011412 < 0.1 Salmon
## 6 NA > 0.1 Salmon
## GENEID baseMean log2FoldChange lfcSE stat pvalue
## 1 ENSG00000223972 16.382284 0.2105319 0.5783491 0.3640221 7.158415e-01
## 2 ENSG00000227232 39.073336 0.2487484 0.3713604 0.6698301 5.029661e-01
## 3 ENSG00000278267 0.000000 NA NA NA NA
## 4 ENSG00000243485 6.640203 0.5409279 0.8752465 0.6180293 5.365560e-01
## 5 ENSG00000284332 0.000000 NA NA NA NA
## 6 ENSG00000237613 202.330477 1.5163697 0.2143040 7.0757887 1.486011e-12
## padj FDR Alignment
## 1 8.876923e-01 > 0.1 STAR
## 2 7.609714e-01 > 0.1 STAR
## 3 NA > 0.1 STAR
## 4 NA > 0.1 STAR
## 5 NA > 0.1 STAR
## 6 1.128834e-10 < 0.1 STAR
## GENEID baseMean log2FoldChange lfcSE stat pvalue
## 1 ENSG00000223972 12.134455 0.42992171 0.6624624 0.64897530 5.163543e-01
## 2 ENSG00000227232 47.652923 -0.02201754 0.3417785 -0.06442051 9.486354e-01
## 3 ENSG00000278267 0.000000 NA NA NA NA
## 4 ENSG00000243485 4.601038 0.83527457 1.0474104 0.79746639 4.251802e-01
## 5 ENSG00000284332 0.000000 NA NA NA NA
## 6 ENSG00000237613 167.762489 1.55049410 0.2207492 7.02378057 2.159439e-12
## padj FDR Alignment
## 1 7.804067e-01 > 0.1 HISAT2
## 2 9.834913e-01 > 0.1 HISAT2
## 3 NA > 0.1 HISAT2
## 4 NA > 0.1 HISAT2
## 5 NA > 0.1 HISAT2
## 6 1.632049e-10 < 0.1 HISAT2# Initialize a data frame storing total lfc results across the aligners
lfc.dataframe <- lfcList[[1]]
for (x in Aligners[2:length(Aligners)]) {
lfc.dataframe <- rbind(lfc.dataframe,
lfcList[[x]])
}
lfc.dataframe$Alignment <- factor(lfc.dataframe$Alignment,
levels=Aligners)Exploring distribution of false discovery rate (FDR)
Black dashed line: FDR = 0.1
# Plot distribution of FDR
ggplot(lfc.dataframe,
aes(x=padj, y=..count.., color=Alignment)) +
geom_density(size=1) +
theme_bw() +
scale_x_log10() +
ggtitle("Distribution of False Discovery Rate (FDR) by Aligner") +
ylab("Count") +
xlim(0.00001, 1) +
geom_vline(xintercept=alpha,
color="black",
size=1, linetype="dashed") + scale_x_continuous(breaks=seq(0, 1, by=0.1))
Presenting distribution of log2FoldChange
Black: total genes (padj =/= NA)
Colored: genes above or below FDR=0.1
valid.lfc.df <- subset(lfc.dataframe, FDR == "< 0.1")
ggplot(valid.lfc.df,
aes(x=log2FoldChange,
y=..count..,
color=Alignment)) +
geom_density(size=1) +
theme_bw() +
geom_vline(xintercept=c(-1, 1),
linetype="dashed", color="black", size=1) +
ggtitle("Distribution of log2FoldChange Values by Aligner (FDR < 0.1)") +
ylab("Count") +
xlim(-10, 10) # Change xlim by datatype
Exploring mean-difference with an MA plot
- Reference: DESeq2 doc “MA-plot”
- x-axis: expression level (baseMean))
- y-axis: fold change (log2FoldChange)
- Red dashed lines: log2FoldChange = -1 and 1
# Set ylim: has to adjusted by users depending on data
yl <- c(-12, 12)
# Set min log2 fold change of interest
mLog <- c(-1, 1)
# Create MA plots by Aligner
ggplot(lfc.dataframe, aes(x=baseMean, y=log2FoldChange, color=FDR)) +
geom_point() +
facet_grid(~Alignment) +
scale_x_log10() +
theme_bw() +
scale_color_manual(values=c("blue", "grey")) +
ggtitle(paste("MA plot")) +
ylim(yl[1], yl[2]) +
theme(strip.text.x=element_text(size=10)) +
geom_hline(yintercept=c(mLog[1], mLog[2]),
linetype="dashed", color="red") 
Exploring expression profiling with normalized count data
- Normalized count matrices are extracted from dds objects and filtered with thresholds set at FDR and log2FoldChange
- The heatmaps display z-scores of the normalized counts
- In this analysis, mLog = 1
- References: Harvard Chan Bioinformatics Core workshop, DESeq2 doc “Heatmap of the count matrix”
# Initialize a list
heatmap.df.List <- lfcList
# Filter genes with FDR < alpha and absolute log2FoldChange > 1
for (x in Aligners) {
# Set a logical vector filtering FDR below alpha
is.fdr.valid <- lfcList[[x]]$FDR == paste("<", alpha)
# Set a logical vector filtering absolute lfc above 1
is.lfc.large <- abs(lfcList[[x]]$log2FoldChange) > mLog[2]
# Extract total normalized counts
norm.counts <- counts(ddsList[[x]], normalized=T)
# Save filtered genes only from the normalized count data
heatmap.df.List[[x]] <- norm.counts[is.fdr.valid & is.lfc.large,]
}
# Explore the cleaned data frames
head(heatmap.df.List[[1]])## Mock-rep1 Mock-rep2 Mock-rep3 SARS-CoV-2-rep1 SARS-CoV-2-rep2
## ENSG00000002746 42.915749 35.23561 27.007451 82.939673 55.73094
## ENSG00000004399 6.391707 12.38008 9.002484 35.964991 26.53854
## ENSG00000005884 61.177769 78.08974 82.147664 330.290732 281.30858
## ENSG00000006327 4.565505 11.42777 15.754346 64.590188 71.65407
## ENSG00000006747 627.300411 690.42755 697.692487 1454.013199 1601.60118
## ENSG00000007001 39.263344 32.37867 34.884624 8.807753 11.94235
## SARS-CoV-2-rep3
## ENSG00000002746 75.29910
## ENSG00000004399 35.20477
## ENSG00000005884 185.80298
## ENSG00000006327 46.93970
## ENSG00000006747 2004.71633
## ENSG00000007001 17.60239head(heatmap.df.List[[2]])## Mock-rep1 Mock-rep2 Mock-rep3 SARS-CoV-2-rep1
## ENSG00000237613 296.941593 294.510344 304.909944 74.09817
## ENSG00000177133 170.718714 143.940747 146.808491 33.01404
## ENSG00000097021 34.506974 35.038208 35.008179 129.85521
## ENSG00000049249 1.816157 14.204679 1.129296 96.10753
## ENSG00000284747 34.506974 40.720080 27.103106 172.40663
## ENSG00000284716 2.724235 5.681872 2.258592 13.93926
## SARS-CoV-2-rep2 SARS-CoV-2-rep3
## ENSG00000237613 129.79225 113.730562
## ENSG00000177133 48.17032 67.650075
## ENSG00000097021 149.86322 86.278357
## ENSG00000049249 104.36903 26.471769
## ENSG00000284747 148.52515 125.495792
## ENSG00000284716 21.40903 7.843487head(heatmap.df.List[[3]])## Mock-rep1 Mock-rep2 Mock-rep3 SARS-CoV-2-rep1 SARS-CoV-2-rep2
## ENSG00000237613 244.74613 244.27485 258.489470 64.00598 108.24497
## ENSG00000177133 163.16409 135.49620 139.875254 32.37084 45.43616
## ENSG00000158286 58.66574 33.39695 40.284073 24.27813 26.72715
## ENSG00000097021 32.08283 32.44275 29.094053 125.80486 134.97212
## ENSG00000049249 0.00000 14.31298 2.238004 93.43402 101.56318
## ENSG00000284747 32.99948 39.12214 22.380041 165.53271 144.32663
## SARS-CoV-2-rep3
## ENSG00000237613 86.81353
## ENSG00000177133 63.40314
## ENSG00000158286 13.65606
## ENSG00000097021 79.98550
## ENSG00000049249 21.45953
## ENSG00000284747 119.97825dim(heatmap.df.List[[1]])## [1] 1087 6dim(heatmap.df.List[[2]])## [1] 1244 6dim(heatmap.df.List[[3]])## [1] 1240 6pheatmap(heatmap.df.List[[3]],
annotation=HeatmapAnno,
scale="row",
show_rownames=F)
# Set a function creating a profiling heatmap
profile.heatmap.fn <- function(df, title) {
pheatmap(df,
annotation=HeatmapAnno,
scale="row",
show_rownames=F,
main=paste("Expression Profiling by", title, "(FDR < 0.1, absolute log2FoldChange > 1)"))
}
# Print the expression heatmaps
profile.heatmap.fn(heatmap.df.List[[Aligners[1]]], Aligners[1])
profile.heatmap.fn(heatmap.df.List[[Aligners[2]]], Aligners[2])
profile.heatmap.fn(heatmap.df.List[[Aligners[3]]], Aligners[3])
NA statistics: zero count genes & outlier genes
When NAs appear in
- log2FoldChange: zero counts in all samples
- padj: too little information
- pval & padj: at least one replicate was an outlier
# Count number of NA genes
type=c("Zero Counts", "Outliers", "Total NA Genes")
# Create a data frame storing number of NA genes by type
NA.genes <- lfc.dataframe %>%
group_by(Alignment) %>%
summarize(zero=sum(is.na(log2FoldChange)),
outlier=sum(is.na(pvalue) & is.na(padj))) %>%
mutate(total=zero + outlier) %>%
gather(Type, Number, -Alignment) %>%
mutate(Type=factor(case_when(Type == "zero" ~ type[1],
Type == "outlier" ~ type[2],
Type == "total" ~ type[3]),
levels=type))
# Plot number of NA genes
ggplot(NA.genes,
aes(x=Type, y=Number, group=Alignment, fill=Alignment, label=Number)) +
geom_bar(stat="identity", position="dodge") +
theme_bw() +
geom_text(position=position_dodge(width=1), vjust=1.5) +
ggtitle("Number of NA Genes") +
ylab("Number of Genes")
Ranking DEGs by alginer
- fdr.rank: ranked by FDR
- lfc.rank: ranked by absolute fold change
- up.lfc.rank: ranked by magnitude of fold change increase
- down.lfc.rank: ranked by manitude of fold change decrease
# Create a new list having DE table with FDR below alpha
fdr.rank <- lfcList
lfc.rank <- lfcList
up.lfc.rank <- lfcList
down.lfc.rank <- lfcList
# Set a sorting genes with FDR below alpha
filter.fdr.fn <- function(df) {as.data.table(df[df$FDR == paste("<", alpha),])}
# Set a function creating a column assigning ranking
Ranking.fn <- function(x) {mutate(x, Rank=1:nrow(x))}
for (x in Aligners) {
rdf <- lfcList[[x]]
fdr.rank[[x]] <- filter.fdr.fn(rdf) %>% arrange(padj) %>% Ranking.fn()
lfc.rank[[x]] <- filter.fdr.fn(rdf) %>% arrange(desc(abs(log2FoldChange))) %>% Ranking.fn()
up.lfc.rank[[x]] <- filter.fdr.fn(rdf) %>% arrange(desc(log2FoldChange)) %>% Ranking.fn()
down.lfc.rank[[x]] <- filter.fdr.fn(rdf) %>% arrange(log2FoldChange) %>% Ranking.fn()
}
# Explore the ranking outputs
head(fdr.rank[[1]])## GENEID baseMean log2FoldChange lfcSE stat pvalue
## 1: ENSG00000261008 1918.3360 -4.546970 0.1938223 -23.45948 1.058028e-121
## 2: ENSG00000169908 2252.9670 -4.043316 0.1751220 -23.08857 6.031531e-118
## 3: ENSG00000118523 3801.4187 -3.486809 0.1825490 -19.10067 2.492575e-81
## 4: ENSG00000136943 11348.9380 -3.121610 0.1746490 -17.87362 1.892955e-71
## 5: ENSG00000181104 2469.3136 -1.949877 0.1160370 -16.80393 2.284141e-63
## 6: ENSG00000109205 356.1611 -4.797026 0.2866798 -16.73305 7.528656e-63
## padj FDR Alignment Rank
## 1: 1.889744e-117 < 0.1 Salmon 1
## 2: 5.386459e-114 < 0.1 Salmon 2
## 3: 1.483996e-77 < 0.1 Salmon 3
## 4: 8.452518e-68 < 0.1 Salmon 4
## 5: 8.159407e-60 < 0.1 Salmon 5
## 6: 2.241155e-59 < 0.1 Salmon 6head(lfc.rank[[1]])## GENEID baseMean log2FoldChange lfcSE stat pvalue
## 1: ENSG00000244398 15.25966 -20.740464 3.9097078 -5.304863 1.127577e-07
## 2: ENSG00000169085 30.97084 -8.403064 1.8213404 -4.613670 3.956199e-06
## 3: ENSG00000283809 30.46514 8.396325 1.2554940 6.687666 2.267577e-11
## 4: ENSG00000198796 51.48336 -8.174275 1.5134337 -5.401145 6.621683e-08
## 5: ENSG00000248167 14.72873 -7.330312 1.3318446 -5.503879 3.715245e-08
## 6: ENSG00000122641 124.98025 -5.844255 0.8233355 -7.098266 1.263314e-12
## padj FDR Alignment Rank
## 1: 3.151745e-06 < 0.1 Salmon 1
## 2: 7.368266e-05 < 0.1 Salmon 2
## 3: 1.281683e-09 < 0.1 Salmon 3
## 4: 1.954874e-06 < 0.1 Salmon 4
## 5: 1.146079e-06 < 0.1 Salmon 5
## 6: 8.579489e-11 < 0.1 Salmon 6head(up.lfc.rank[[1]])## GENEID baseMean log2FoldChange lfcSE stat pvalue
## 1: ENSG00000283809 30.465135 8.396325 1.2554940 6.687666 2.267577e-11
## 2: ENSG00000271723 19.808302 4.756011 0.9297970 5.115107 3.135631e-07
## 3: ENSG00000265168 17.736668 4.596849 0.9245614 4.971924 6.629158e-07
## 4: ENSG00000279750 9.865767 3.825804 1.0778120 3.549603 3.858130e-04
## 5: ENSG00000251177 8.980151 3.764548 1.0508629 3.582340 3.405298e-04
## 6: ENSG00000105641 14.315585 3.174151 0.7758214 4.091342 4.288844e-05
## padj FDR Alignment Rank
## 1: 1.281683e-09 < 0.1 Salmon 1
## 2: 7.982744e-06 < 0.1 Salmon 2
## 3: 1.535712e-05 < 0.1 Salmon 3
## 4: 3.862671e-03 < 0.1 Salmon 4
## 5: 3.475137e-03 < 0.1 Salmon 5
## 6: 6.026990e-04 < 0.1 Salmon 6head(down.lfc.rank[[1]])## GENEID baseMean log2FoldChange lfcSE stat pvalue
## 1: ENSG00000244398 15.25966 -20.740464 3.9097078 -5.304863 1.127577e-07
## 2: ENSG00000169085 30.97084 -8.403064 1.8213404 -4.613670 3.956199e-06
## 3: ENSG00000198796 51.48336 -8.174275 1.5134337 -5.401145 6.621683e-08
## 4: ENSG00000248167 14.72873 -7.330312 1.3318446 -5.503879 3.715245e-08
## 5: ENSG00000122641 124.98025 -5.844255 0.8233355 -7.098266 1.263314e-12
## 6: ENSG00000187689 59.70984 -5.675160 1.1955478 -4.746912 2.065462e-06
## padj FDR Alignment Rank
## 1: 3.151745e-06 < 0.1 Salmon 1
## 2: 7.368266e-05 < 0.1 Salmon 2
## 3: 1.954874e-06 < 0.1 Salmon 3
## 4: 1.146079e-06 < 0.1 Salmon 4
## 5: 8.579489e-11 < 0.1 Salmon 5
## 6: 4.094474e-05 < 0.1 Salmon 6Calculating rank difference between STAR and HISAT2
Salmon was excluded due to extremely small number of genes found
# Set a function rebuilding DE tables with gene ranks
rankdiff.fn <- function(List){
# Select columns and join the data frames by gene
full_join(List[[Aligners[2]]][, .(GENEID, Alignment, Rank, baseMean)],
List[[Aligners[3]]][, .(GENEID, Alignment, Rank, baseMean)],
by="GENEID") %>%
# Add columns assining gene expression levels, rank differences, and mean ranks
mutate(logMeanExpression=log(baseMean.x+baseMean.y/2),
RankDiff=Rank.x-Rank.y,
MeanRank=(Rank.x+Rank.y)/2)
}
# Calculate rank difference by ranking type
fdr.rankdiff <- rankdiff.fn(fdr.rank)
lfc.rankdiff <- rankdiff.fn(lfc.rank)
up.lfc.rankdiff <- rankdiff.fn(up.lfc.rank)
down.lfc.rankdiff <- rankdiff.fn(down.lfc.rank)
# Explore the calculated rank differences
head(fdr.rankdiff)## GENEID Alignment.x Rank.x baseMean.x Alignment.y Rank.y baseMean.y
## 1: ENSG00000261008 STAR 1 2036.3102 HISAT2 1 1977.5615
## 2: ENSG00000169908 STAR 2 2336.0090 HISAT2 2 2262.3133
## 3: ENSG00000118523 STAR 3 3921.7264 HISAT2 3 3820.8621
## 4: ENSG00000136943 STAR 4 11817.9679 HISAT2 4 11418.8145
## 5: ENSG00000109205 STAR 5 382.4466 HISAT2 5 374.5982
## 6: ENSG00000181104 STAR 6 2551.1683 HISAT2 6 2478.0294
## logMeanExpression RankDiff MeanRank
## 1: 8.014696 0 1
## 2: 8.151093 0 2
## 3: 8.671142 0 3
## 4: 9.771519 0 4
## 5: 6.345190 0 5
## 6: 8.240170 0 6head(lfc.rankdiff)## GENEID Alignment.x Rank.x baseMean.x Alignment.y Rank.y baseMean.y
## 1: ENSG00000259895 STAR 1 200.612721 <NA> NA NA
## 2: ENSG00000138675 STAR 2 41.175215 HISAT2 3 39.989193
## 3: ENSG00000178776 STAR 3 8.745942 HISAT2 2 8.182662
## 4: ENSG00000198796 STAR 4 37.920036 HISAT2 1 37.316116
## 5: ENSG00000187689 STAR 5 64.800783 HISAT2 5 62.192291
## 6: ENSG00000122641 STAR 6 128.103619 HISAT2 4 124.739803
## logMeanExpression RankDiff MeanRank
## 1: NA NA NA
## 2: 4.113654 -1 2.5
## 3: 2.552353 1 2.5
## 4: 4.035622 3 2.5
## 5: 4.563274 0 5.0
## 6: 5.249513 2 5.0head(up.lfc.rankdiff)## GENEID Alignment.x Rank.x baseMean.x Alignment.y Rank.y baseMean.y
## 1: ENSG00000259895 STAR 1 200.612721 <NA> NA NA
## 2: ENSG00000103184 STAR 2 8.588905 HISAT2 8 7.399501
## 3: ENSG00000105641 STAR 3 14.872126 HISAT2 1 13.620723
## 4: ENSG00000256654 STAR 4 15.529629 HISAT2 7 14.733459
## 5: ENSG00000251177 STAR 5 9.858227 HISAT2 5 9.732927
## 6: ENSG00000144649 STAR 6 7.421817 HISAT2 2 6.526874
## logMeanExpression RankDiff MeanRank
## 1: NA NA NA
## 2: 2.508677 -6 5.0
## 3: 3.076505 2 2.0
## 4: 3.130978 -3 5.5
## 5: 2.689526 0 5.0
## 6: 2.368865 4 4.0head(down.lfc.rankdiff)## GENEID Alignment.x Rank.x baseMean.x Alignment.y Rank.y baseMean.y
## 1: ENSG00000138675 STAR 1 41.175215 HISAT2 3 39.989193
## 2: ENSG00000178776 STAR 2 8.745942 HISAT2 2 8.182662
## 3: ENSG00000198796 STAR 3 37.920036 HISAT2 1 37.316116
## 4: ENSG00000187689 STAR 4 64.800783 HISAT2 5 62.192291
## 5: ENSG00000122641 STAR 5 128.103619 HISAT2 4 124.739803
## 6: ENSG00000137673 STAR 6 143.085583 HISAT2 7 137.747910
## logMeanExpression RankDiff MeanRank
## 1: 4.113654 -2 2.0
## 2: 2.552353 0 2.0
## 3: 4.035622 2 2.0
## 4: 4.563274 -1 4.5
## 5: 5.249513 1 4.5
## 6: 5.356395 -1 6.5Visualizing DEG rankings
- x-axis: genes aligned by STAR
- y-axis: genes aligned by HISAT2
- Black diagonal lines: ranking in STAR = HISAT2
- Dot color: gene expression level (log-baseMean)
- 409 genes were missing in the plots
# Set a function plotting DEG rankings
ranking.plot.fn <- function(df, rankedby) {
ggplot(df,
aes(x=Rank.x, y=Rank.y, color=logMeanExpression)) + geom_point(alpha=0.5) + theme_bw() + theme(strip.text.x=element_text(size=10)) + xlab("Ranking in STAR") + ylab("Ranking in HISAT2") + geom_abline(slope=1, color="black", size=0.5) + ggtitle(paste("Gene Ranking by", rankedby)) + scale_color_gradient(low="blue", high="red")
}
# Plot rankings by ranking type
ranking.plot.fn(fdr.rankdiff, "FDR")
ranking.plot.fn(lfc.rankdiff, "Log2FoldChange")
ranking.plot.fn(up.lfc.rankdiff, "Log2FoldChange (Increase)")
ranking.plot.fn(down.lfc.rankdiff, "Log2FoldChange (Decrease)")
Visualizing DEG ranks II: Rank difference
- x-axis: expression level (log-baseMean)
- y-axis: rank difference (STAR ranking - HISAT2 ranking)
- Black horizontal lines: ranking in STAR = HISAT2
- Dot color: average ranking
- 409 genes were missing in the plots
# Set a function plotting DEG rank difference
rankdiff.plot.fn <- function(df, rankedby, ylim) {
ggplot(df, aes(x=logMeanExpression, y=RankDiff, color=MeanRank)) +
geom_point(alpha=0.5) +
theme_bw() +
ylab("Rank Difference (STAR - HISAT2)") +
ggtitle(paste("Rank Difference (STAR - HISAT2)\nin", rankedby)) +
geom_hline(yintercept=0, color="black", size=0.5) + scale_color_gradient(low="blue", high="red") +
ylim(-ylim, ylim)
}
# Display the plots in the same y-scale
rankdiff.plot.fn(fdr.rankdiff, "FDR", 2500) 
rankdiff.plot.fn(lfc.rankdiff, "Log2FoldChange", 2500)
rankdiff.plot.fn(up.lfc.rankdiff, "Log2FoldChange (Increase)", 2500)
rankdiff.plot.fn(down.lfc.rankdiff, "Log2FoldChange (Decrease)", 2500)
# Display the plots in free y-scale
rankdiff.plot.fn(lfc.rankdiff, "Log2FoldChange (free y-scale)", 1500)
rankdiff.plot.fn(up.lfc.rankdiff, "Log2FoldChange (Increase, free y-scale)", 750)
rankdiff.plot.fn(down.lfc.rankdiff, "Log2FoldChange (Decrease, free y-scale)", 600)
Distribution of rank difference
# Create a list storing rankdiff data frames
rankList <- list(fdr.rankdiff,
lfc.rankdiff,
up.lfc.rankdiff,
down.lfc.rankdiff)
# Assine result column as a factor with levels
rankdiff.levels <- c("FDR",
"log2FoldChange",
"log2FoldChange.Increase",
"log2FoldChange.Decrease")
# Name the list
names(rankList) <- rankdiff.levels
# Create a new data frame storing rank difference by result type
rankdiff.dist <- data.frame(FDR=abs(rankList[[1]]$RankDiff),
log2FoldChange=abs(rankList[[2]]$RankDiff),
log2FoldChange.Increase=abs(rankList[[3]]$RankDiff),
log2FoldChange.Decrease=abs(rankList[[4]]$RankDiff)) %>% gather(Result, RankDiff)
rankdiff.dist$Result <- factor(rankdiff.dist$Result, levels=rankdiff.levels)
# Plot distribution of absolute rank difference
ggplot(rankdiff.dist,
aes(x=Result, y=RankDiff, color=Result)) +
geom_jitter(alpha=0.5) +
geom_boxplot(alpha=0.5, fill="grey", color="black") +
theme_bw() +
theme(axis.text.x=element_text(angle=45, hjust=1)) +
ggtitle("Distribution of Absolute Rank Difference without Shrinkage \n(STAR - HISAT2)") +
ylab("Absolute Rank Difference (STAR - HISAT2)") 
Relationship between rank difference and number of transcript versions
- y-axis: abs(TPM-Count inputs)
- x-axis: number of transcripts (number of transcript id / gene id)
- dot color: mean rank
# Create a data frame storing the number of transcripts by gene id
AnnoDb.ntrans <- AnnoDb %>%
group_by(GENEID) %>%
summarize(num.trans=n_distinct(TXID))
# Set a function adding the number of transcripts by gene id
add.ntrans.fn <- function(df) {
left_join(df, AnnoDb.ntrans, by="GENEID")}
# Add a column indicating the number of transcripts by gene id to every rankdiff data frame
for (x in rankdiff.levels) {
rankList[[x]] <- add.ntrans.fn(rankList[[x]])
}
# Explore the edited data frames
summary(rankList)## Length Class Mode
## FDR 11 data.table list
## log2FoldChange 11 data.table list
## log2FoldChange.Increase 11 data.table list
## log2FoldChange.Decrease 11 data.table listhead(rankList[[1]])## GENEID Alignment.x Rank.x baseMean.x Alignment.y Rank.y baseMean.y
## 1: ENSG00000261008 STAR 1 2036.3102 HISAT2 1 1977.5615
## 2: ENSG00000169908 STAR 2 2336.0090 HISAT2 2 2262.3133
## 3: ENSG00000118523 STAR 3 3921.7264 HISAT2 3 3820.8621
## 4: ENSG00000136943 STAR 4 11817.9679 HISAT2 4 11418.8145
## 5: ENSG00000109205 STAR 5 382.4466 HISAT2 5 374.5982
## 6: ENSG00000181104 STAR 6 2551.1683 HISAT2 6 2478.0294
## logMeanExpression RankDiff MeanRank num.trans
## 1: 8.014696 0 1 18
## 2: 8.151093 0 2 4
## 3: 8.671142 0 3 1
## 4: 9.771519 0 4 3
## 5: 6.345190 0 5 5
## 6: 8.240170 0 6 2head(rankList[[2]])## GENEID Alignment.x Rank.x baseMean.x Alignment.y Rank.y baseMean.y
## 1: ENSG00000259895 STAR 1 200.612721 <NA> NA NA
## 2: ENSG00000138675 STAR 2 41.175215 HISAT2 3 39.989193
## 3: ENSG00000178776 STAR 3 8.745942 HISAT2 2 8.182662
## 4: ENSG00000198796 STAR 4 37.920036 HISAT2 1 37.316116
## 5: ENSG00000187689 STAR 5 64.800783 HISAT2 5 62.192291
## 6: ENSG00000122641 STAR 6 128.103619 HISAT2 4 124.739803
## logMeanExpression RankDiff MeanRank num.trans
## 1: NA NA NA 1
## 2: 4.113654 -1 2.5 5
## 3: 2.552353 1 2.5 3
## 4: 4.035622 3 2.5 6
## 5: 4.563274 0 5.0 2
## 6: 5.249513 2 5.0 5head(rankList[[3]])## GENEID Alignment.x Rank.x baseMean.x Alignment.y Rank.y baseMean.y
## 1: ENSG00000259895 STAR 1 200.612721 <NA> NA NA
## 2: ENSG00000103184 STAR 2 8.588905 HISAT2 8 7.399501
## 3: ENSG00000105641 STAR 3 14.872126 HISAT2 1 13.620723
## 4: ENSG00000256654 STAR 4 15.529629 HISAT2 7 14.733459
## 5: ENSG00000251177 STAR 5 9.858227 HISAT2 5 9.732927
## 6: ENSG00000144649 STAR 6 7.421817 HISAT2 2 6.526874
## logMeanExpression RankDiff MeanRank num.trans
## 1: NA NA NA 1
## 2: 2.508677 -6 5.0 2
## 3: 3.076505 2 2.0 2
## 4: 3.130978 -3 5.5 20
## 5: 2.689526 0 5.0 1
## 6: 2.368865 4 4.0 5head(rankList[[4]])## GENEID Alignment.x Rank.x baseMean.x Alignment.y Rank.y baseMean.y
## 1: ENSG00000138675 STAR 1 41.175215 HISAT2 3 39.989193
## 2: ENSG00000178776 STAR 2 8.745942 HISAT2 2 8.182662
## 3: ENSG00000198796 STAR 3 37.920036 HISAT2 1 37.316116
## 4: ENSG00000187689 STAR 4 64.800783 HISAT2 5 62.192291
## 5: ENSG00000122641 STAR 5 128.103619 HISAT2 4 124.739803
## 6: ENSG00000137673 STAR 6 143.085583 HISAT2 7 137.747910
## logMeanExpression RankDiff MeanRank num.trans
## 1: 4.113654 -2 2.0 5
## 2: 2.552353 0 2.0 3
## 3: 4.035622 2 2.0 6
## 4: 4.563274 -1 4.5 2
## 5: 5.249513 1 4.5 5
## 6: 5.356395 -1 6.5 3# Set a function plotting rank difference vs number of transcripts
rank.ntrans.plot.fn <- function(df, title) {
ggplot(df, aes(x=num.trans, y=abs(RankDiff), color=MeanRank)) +
geom_jitter(alpha=0.5) +
theme_bw() +
ggtitle(paste("Rank Difference vs Number of Alternative Transcripts \nin", title)) +
xlab("Number of Alternative Transcripts") +
ylab("Absolute Rank Difference \n(STAR - HISAT2)") + scale_color_gradient(low="blue", high="red")
}
# Print the plots
rank.ntrans.plot.fn(rankList[[1]], "FDR")
rank.ntrans.plot.fn(rankList[[2]], "log2FoldChange")
rank.ntrans.plot.fn(rankList[[3]], "log2FoldChange (Increase)")
rank.ntrans.plot.fn(rankList[[4]], "log2FoldChange (Decrease)")
Finding genes having large difference in rankings
# Initialize a list storing rankdiff genes
large.rankdiff <- rankList
# Assign a vector storing minimum (thresholds) rankdiff for filtering large rankdiff genes
rankdiff.thr <- c(10, 10, 10, 10)
names(rankdiff.thr) <- rankdiff.levels
for (x in rankdiff.levels) {
# Filter out observations below the rankdiff thresholds
large.rankdiff[[x]] <- subset(rankList[[x]],
abs(RankDiff) > rankdiff.thr[x]) %>%
# Arrange by descending order of RankDiff
arrange(desc(abs(RankDiff)))
}
# Explore the filtered genes
summary(large.rankdiff)## Length Class Mode
## FDR 11 data.table list
## log2FoldChange 11 data.table list
## log2FoldChange.Increase 11 data.table list
## log2FoldChange.Decrease 11 data.table listdim(large.rankdiff[[rankdiff.levels[1]]])## [1] 3631 11dim(large.rankdiff[[rankdiff.levels[2]]])## [1] 3637 11dim(large.rankdiff[[rankdiff.levels[3]]])## [1] 3613 11dim(large.rankdiff[[rankdiff.levels[4]]])## [1] 3733 11head(large.rankdiff[[rankdiff.levels[1]]])## GENEID Alignment.x Rank.x baseMean.x Alignment.y Rank.y baseMean.y
## 1: ENSG00000151835 STAR 1040 2773.89592 HISAT2 3329 2718.41143
## 2: ENSG00000215630 STAR 1995 202.68947 HISAT2 4203 112.14843
## 3: ENSG00000168333 STAR 2966 93.73031 HISAT2 822 89.42132
## 4: ENSG00000198535 STAR 2684 32.21448 HISAT2 571 31.83286
## 5: ENSG00000280048 STAR 753 260.63943 HISAT2 2780 252.81904
## 6: ENSG00000134184 STAR 2255 199.52666 HISAT2 4253 148.34171
## logMeanExpression RankDiff MeanRank num.trans
## 1: 8.326783 -2289 2184.5 5
## 2: 5.555915 -2208 3099.0 1
## 3: 4.930444 2144 1894.0 4
## 4: 3.873925 2113 1627.5 1
## 5: 5.958551 -2027 1766.5 1
## 6: 5.612024 -1998 3254.0 8head(large.rankdiff[[rankdiff.levels[2]]])## GENEID Alignment.x Rank.x baseMean.x Alignment.y Rank.y baseMean.y
## 1: ENSG00000282508 STAR 2972 195.96247 HISAT2 1845 73.61686
## 2: ENSG00000244722 STAR 2126 94.85277 HISAT2 1162 65.13034
## 3: ENSG00000241549 STAR 2808 296.10784 HISAT2 1866 147.00192
## 4: ENSG00000196302 STAR 2440 180.23862 HISAT2 1521 92.31848
## 5: ENSG00000213569 STAR 2021 166.51235 HISAT2 2869 186.29704
## 6: ENSG00000267544 STAR 2666 217.56722 HISAT2 1879 194.83397
## logMeanExpression RankDiff MeanRank num.trans
## 1: 5.450055 1127 2408.5 39
## 2: 4.847473 964 1644.0 1
## 3: 5.912445 942 2337.0 3
## 4: 5.422294 919 1980.5 1
## 5: 5.559376 -848 2445.0 1
## 6: 5.752523 787 2272.5 1head(large.rankdiff[[rankdiff.levels[3]]])## GENEID Alignment.x Rank.x baseMean.x Alignment.y Rank.y baseMean.y
## 1: ENSG00000282508 STAR 1386 195.96247 HISAT2 830 73.61686
## 2: ENSG00000213569 STAR 3324 166.51235 HISAT2 2806 186.29704
## 3: ENSG00000241549 STAR 1298 296.10784 HISAT2 842 147.00192
## 4: ENSG00000196302 STAR 1106 180.23862 HISAT2 658 92.31848
## 5: ENSG00000244722 STAR 3270 94.85277 HISAT2 3623 65.13034
## 6: ENSG00000172345 STAR 4322 15.63896 HISAT2 4027 28.17195
## logMeanExpression RankDiff MeanRank num.trans
## 1: 5.450055 556 1108.0 39
## 2: 5.559376 518 3065.0 1
## 3: 5.912445 456 1070.0 3
## 4: 5.422294 448 882.0 1
## 5: 4.847473 -353 3446.5 1
## 6: 3.391986 295 4174.5 6head(large.rankdiff[[rankdiff.levels[4]]])## GENEID Alignment.x Rank.x baseMean.x Alignment.y Rank.y baseMean.y
## 1: ENSG00000244722 STAR 1179 94.85277 HISAT2 697 65.13034
## 2: ENSG00000282508 STAR 3063 195.96247 HISAT2 3490 73.61686
## 3: ENSG00000267544 STAR 1448 217.56722 HISAT2 1033 194.83397
## 4: ENSG00000230847 STAR 3729 633.44107 HISAT2 3334 551.50468
## 5: ENSG00000213569 STAR 1125 166.51235 HISAT2 1514 186.29704
## 6: ENSG00000232037 STAR 2855 1717.75938 HISAT2 2470 1441.56074
## logMeanExpression RankDiff MeanRank num.trans
## 1: 4.847473 482 938.0 1
## 2: 5.450055 -427 3276.5 39
## 3: 5.752523 415 1240.5 1
## 4: 6.812558 395 3531.5 1
## 5: 5.559376 -389 1319.5 1
## 6: 7.799155 385 2662.5 1# Set a function saving rankdiff.csv files
saving.fn <- function(input.df, data.type) {
filename <- paste0(DB, "_DE_", data.type, "_RankDiff.csv")
return(write.csv(input.df, filename))
}
# Save the filtered and arranged data frames as csv files
saving.fn(large.rankdiff[[rankdiff.levels[1]]], "FDR")
saving.fn(large.rankdiff[[rankdiff.levels[2]]], "LFC")
saving.fn(large.rankdiff[[rankdiff.levels[3]]], "LFC_Increase")
saving.fn(large.rankdiff[[rankdiff.levels[4]]], "LFC_Decrease")Summarizing up/down DEGs with an upset plot
red bar: aligner
blue bar: directionality of gene expression change
# Generate a data frame storing upset input variables
upset.dataframe <- subset(lfc.dataframe, !is.na(padj)) %>%
mutate(Up=ifelse(FDR == paste("<", alpha) & log2FoldChange > 0, GENEID, ""),
Down=ifelse(FDR == paste("<", alpha) & log2FoldChange < 0, GENEID, ""),
Unchanged=ifelse(FDR == paste(">", alpha), GENEID, ""),
Salmon=ifelse(Alignment == "Salmon", GENEID, ""),
STAR=ifelse(Alignment == "STAR", GENEID, ""),
HISAT2=ifelse(Alignment == "HISAT2", GENEID, ""))
# Generate a list input
upset.input <- list(Up=upset.dataframe$Up,
Down=upset.dataframe$Down,
Unchanged=upset.dataframe$Unchanged,
Salmon=upset.dataframe$Salmon,
STAR=upset.dataframe$STAR,
HISAT2=upset.dataframe$HISAT2)
# Create an upset plot
upset(fromList(upset.input),
sets=names(upset.input), # What group to display
sets.x.label="Number of Genes per Group",
order.by="freq",
point.size=3,
sets.bar.color=c("red", "red", "blue", "red", "blue", "blue"),
text.scale = 1.5, number.angles=30)
Session Info
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /home/mira/miniconda3/envs/r/lib/libblas.so.3.8.0
## LAPACK: /home/mira/miniconda3/envs/r/lib/liblapack.so.3.8.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] ensembldb_2.14.0 AnnotationFilter_1.14.0
## [3] GenomicFeatures_1.42.0 AnnotationDbi_1.52.0
## [5] UpSetR_1.4.0 DESeq2_1.30.1
## [7] SummarizedExperiment_1.20.0 Biobase_2.50.0
## [9] MatrixGenerics_1.2.0 matrixStats_0.58.0
## [11] GenomicRanges_1.42.0 GenomeInfoDb_1.26.0
## [13] IRanges_2.24.0 S4Vectors_0.28.0
## [15] Rsubread_2.4.0 tximport_1.18.0
## [17] AnnotationHub_2.22.0 BiocFileCache_1.14.0
## [19] dbplyr_2.1.0 BiocGenerics_0.36.0
## [21] pheatmap_1.0.12 rmarkdown_2.7
## [23] forcats_0.5.1 stringr_1.4.0
## [25] dplyr_1.0.4 purrr_0.3.4
## [27] readr_1.4.0 tidyr_1.1.2
## [29] tibble_3.0.6 ggplot2_3.3.2
## [31] tidyverse_1.3.0 data.table_1.13.4
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-0 ellipsis_0.3.1
## [3] XVector_0.30.0 fs_1.5.0
## [5] rstudioapi_0.13 farver_2.0.3
## [7] bit64_4.0.5 interactiveDisplayBase_1.28.0
## [9] fansi_0.4.2 lubridate_1.7.10
## [11] xml2_1.3.2 splines_4.0.3
## [13] cachem_1.0.4 geneplotter_1.68.0
## [15] knitr_1.31 jsonlite_1.7.2
## [17] Rsamtools_2.6.0 broom_0.7.5
## [19] annotate_1.68.0 shiny_1.6.0
## [21] BiocManager_1.30.10 compiler_4.0.3
## [23] httr_1.4.2 backports_1.2.1
## [25] lazyeval_0.2.2 assertthat_0.2.1
## [27] Matrix_1.3-2 fastmap_1.1.0
## [29] cli_2.3.1 later_1.1.0.1
## [31] prettyunits_1.1.1 htmltools_0.5.1.1
## [33] tools_4.0.3 gtable_0.3.0
## [35] glue_1.4.2 GenomeInfoDbData_1.2.4
## [37] rappdirs_0.3.3 Rcpp_1.0.6
## [39] cellranger_1.1.0 jquerylib_0.1.3
## [41] Biostrings_2.58.0 vctrs_0.3.6
## [43] rtracklayer_1.50.0 xfun_0.20
## [45] ps_1.5.0 rvest_0.3.6
## [47] mime_0.10 lifecycle_1.0.0
## [49] XML_3.99-0.5 zlibbioc_1.36.0
## [51] scales_1.1.1 ProtGenerics_1.22.0
## [53] hms_1.0.0 promises_1.2.0.1
## [55] RColorBrewer_1.1-2 yaml_2.2.1
## [57] curl_4.3 gridExtra_2.3
## [59] memoise_2.0.0 sass_0.3.1
## [61] biomaRt_2.46.0 stringi_1.5.3
## [63] RSQLite_2.2.3 highr_0.8
## [65] BiocVersion_3.12.0 genefilter_1.72.1
## [67] BiocParallel_1.24.0 rlang_0.4.10
## [69] pkgconfig_2.0.3 bitops_1.0-6
## [71] evaluate_0.14 lattice_0.20-41
## [73] labeling_0.4.2 GenomicAlignments_1.26.0
## [75] bit_4.0.4 tidyselect_1.1.0
## [77] plyr_1.8.6 magrittr_2.0.1
## [79] R6_2.5.0 generics_0.1.0
## [81] DelayedArray_0.16.0 DBI_1.1.1
## [83] pillar_1.5.0 haven_2.3.1
## [85] withr_2.4.1 survival_3.2-7
## [87] RCurl_1.98-1.2 modelr_0.1.8
## [89] crayon_1.4.1 utf8_1.1.4
## [91] progress_1.2.2 locfit_1.5-9.4
## [93] grid_4.0.3 readxl_1.3.1
## [95] blob_1.2.1 reprex_1.0.0
## [97] digest_0.6.27 xtable_1.8-4
## [99] httpuv_1.5.5 openssl_1.4.3
## [101] munsell_0.5.0 bslib_0.2.4
## [103] askpass_1.1